■a 


'|i  'WIj  mS^  '5 _ *)fk  4* 

;ts  ,,%■■  v.  ■5:  ;m  !•»  5*1  Iki  |ii  8?b  is  11  ir  ■  it  i»  ic  it  !», 


M*  ifi  M‘j*|j'  iT|'r|T 
■7 

esei  "rv\"i2  .ff 


ii 


4» 


Hi 

i 


iiiii  tl{  llliiii!!|i 


n|^e^  ^  •eoMe 


i 

I 


i| 

Ills  III 


1  I 
' 


).u  !i!!<  Itli  liiHui!  ii 


II 


*1 

!? 

I 


f 


if 

i! 

I! 

?l 

}l 


J«lf  'Si'SrYk  8« 


tl«  •ftvVS/V 


secufil^v  CU  ASStrtCATlON  or  This  page  fWh0n  emfi04) 


REPORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


RCC«P|CII|T*|  catalog  NUMOCA 


A.  Title  (and  SuAlirf*; 


Numerical  Implementation  of  the  Cohesive 
Soil  Bounding  Surface  Plasticity  Model 
(Volume  I) 


7  AUTMOAfe; 

Leonard  R.  Herrmann,  Prof.;  Jay  S. 
DeNatale,  Asst.  Prof.;  Yannis  F.  Oafalia: 


%.  TVAC  OP  ACPOAT  A  AEAIOO  CQVCACO 

Not  Final 

Apr  1982  -  Nov  1982 


•  PCAPOAMIMG  OAC.  ACPOAT  MUMOCA 


OMTAACT  OA  CAANT  NUMOCATsi 


N62474-82-C-8276 


*  PCaFORUINC  - - - 

University  of  California 
Davis,  CA 


1  •  controlling  office  name  ano  aooaess 

Naval  Civil  Engineering  Laboratory 
Port  Hueneme,  CA  93043 


T4  MONITOAinC  agency  name  a  AOOACSVif  dUtatant  Imm  C^ntt^Hirng  Otticd) 

Naval  Facilities  Engineering  Command 
200  Stovall  Street 
Alexandria,  VA  22332 


lA  OlSTAlAuTlON  statement  Mapott} 


Approved  for  public  release;  distribution  unlimited 


17  OtSTAiOuTiON  statement  fpt  tha  pbatept*  enraretf  M  RfaeA  7Q.  tt  tram  AepArf) 


i2.  AEPOAT  date 

February  1983 


tS  NUMAEA  OP  PAGES 

114 


II.  lECUAlTY  CLASS  (ot  rMe  r«port> 


Unclassified 


19  NET  WO  AOS  (Cenifnwe  on  reverie  erWe  lliteeeaaery  enP  IPamifr  Af  MacA  mewAer} 


Constitutive  law,  cohesive  soil,  finite  element,  numerical 
implementation 


20^  I 

^The  results  of  a  study  of  various  numerical  means  for  \ 

implementing  the  bounding  surface  plasticity  model  for  cohesive 
soils  is  presented.  A  comparison  is  made  of  the  computational 
efficiency,  robustness, and  ease  of  implementation  of  the  com¬ 
pleting  methods.  The  comparisons  are  made  on  the  basis  of 
solut^ns  to  three  representative  geotechnical  engineering 

00  ,  1473  eo,T,oN  OF .  NOV .. .. oMOLtre  Unclassified 

»«eu"ITT  CtNMiFICATIOM  OF  TNI*  NAC«  !>•••  tnlntd) 


SCCunirv  Ct  a&&*^«Cation  o**  this  PA(ic  riHiM  Dm* 


>  A  very  brief  description  of  a  computer-aided  automated 
calibration  procedure  is  given  along  with  instructions  for  usin 
the  resulting  computer  code.  An  example  calibration  is  report* 
for  a  particular  cohesive  soil. 


.  ! 


. . -.J 


DD  I  jan'ti  1473  COITION  or  I  MOV  ft  If  OOfOLCTC 


Unclassified _ 

$CCuK4f Y  CWMtl^tCATiOii  Of  TMtf  P*OC  Of 


TABLE  OF  CONTENTS 


ABSTRACT  I 

1.  INTRODUCTION  2 

2.  THE  AUTOMATED  CALIBRATION  CODE  3 

2.1  Purpose  and  Capabilities  3 

2.2  The  Calibration  Data  Base  3 

2.3  Practical  Considerations  8 

2.4  Example  9 

2.5  Cost  10 

3.  NUMERICAL  IMPLEMENTATION  10 

3.1  Incrementalization  Of  Rate  Equations  10 

3.2  Calculation  Of  Pore  Water  Pressure  14 

3.3  General  Discussion  Of  Finite  Element  Implementation 

Of  Plasticity  Models  17 

3.4  Selection  Of  Methods  For  Comparison  17 

3.5  Theory  20 

3.5.1  Solution  Methods  20 

3.5.2  Reduction  Of  The  Number  Of  Equation 

Triangularizations  24 

3.5.3  Convergence  Acceleration  27 

3.5.4  Accounting  For  The  Inconsistency  In  The  Incremental 

Stress-Strain  Relation  28 

3.6  Finite  Element  Implementation  29 

4.  NUMERICAL  STUDY  30 

4.1  Scope  Of  Study  30 

4.2  Selection  OI  Representative  Problems  30 

4.3  Comparison  Results  32 

5.  CONCLUSIONS  35 

6.  RECOMMENDED  FUTURE  WORK  36 

REFERENCES  37 

TABLES  40 

FIGURES  42 

APPENDIX  A:  NON-HOMOGENEOUS  MODEL  AND  IN  SITU 

TEST  RESULTS  52 

APPENDIX  B:  USER'S  MANUAL  FOR  NTD  55 


1.  INTRODUCTION 

A  preliminary  task  of  the  project  involved  the  continued  evaluation  of  certain 
features  of  the  model  and  the  reporting  of  any  additional  theoretical  developments. 

Under  other  sponsorship,  Professor  Dafalias  continued  to  study  the  predictive 
capabilities  of  the  model  for  cyclic  soil  behavior.  Preliminary  results  (Dafalias  et  al 
1982)  indicate  the  need  for  a  mechanism  to  model  the  degradation  of  properties  observed 
under  long-term  cyclic  loading  conditions.  A  possible  means  for  modeling  this 
phenomenon  was  developed  and  some  preliminary  testing  was  conducted.  However,  a 
complete  validation  of  this  modeling  mechanism,  and  the  determination  of  its  usefulness, 
will  require  additional  study. 

Continued  compeu-isons  of  model  predictions  and  experimental  results  indicate 
the  desirability  of  using  a  projection  point  other  than  the  origin  for  the  mapping  rule. 
Accordingly,  this  feature  has  been  included  in  the  current  version  of  the  model.  The 
interaction  of  this  feature  with  other  characteristics  of  the  model  is  discussed  in  detail 
by  DeNatale  (1982). 

Numerical  evaluations  of  the  model  made  during  the  course  of  the  current  study 
revealed  some  minor  numerical  problems  associated  with  the  first  solution  step  away 
from  a  pure  hydrostatic  stress  state  for  a  normally  consolidated  soil.  Professor  Dafalias 
has  proposed  a  slight  modification  of  the  model  to  remedy  this  problem.  When  this 
suggested  change  has  been  fully  evaluated,  and  assuming  that  it  does  rectify  the 
problem,  it  will  be  incorporated  into  the  model. 

The  development  of  a  computer  aided,  automated  calibration  scheme  has  been 
completed  and  used  to  perform  a  very  extensive  study  of  the  importance  and  roles  of 
the  several  material  parameters  which  describe  the  model.  A  brief  description  of  this 
study  is  given  in  the  following  section,  and  volume  n  of  this  report  contains  a  user's 
manual  and  a  listing  for  the  calibration  program.  Before  the  end  of  the  year,  the 
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complete  results  of  the  study  will  be  available  in  DeNatale's  thesis  (1982);  a  copy  will 
be  forwarded  to  NCEL  at  that  time. 

The  main  thrust  of  the  current  work  has  been  the  study  of  various  schemes  for 
numerically  implementing  the  model  in  finite  element  analyses  and  the  study  of  the 
numerical  characteristics  of  the  model.  Particular  attention  was  paid  to  the  relative 
ease  of  implementation,  and  the  economy  and  robustness  of  the  competing  schemes. 
In  previous  work  (Herrmann  et  al  1981),  the  rate  equations  for  the  bounding  surface 
model  were  cast  in  incremental  form,  and  a  subroutine  was  prepared  to  evaluate  them. 
In  Section  3,  a  very  brief  review  of  this  step  will  be  given,  prior  to  the  reporting  on 
the  main  component  of  the  research. 

The  final  phase  of  the  overall  project  will  involve  the  verification  of  the 
predictive  capabilities  of  finite  element  analyses  which  utilize  the  model.  With  this 
end  in  mind,  a  preliminary  literature  survey  was  conducted  to  determine  the  availability 
of  laboratory  and  field  experimental  results.  The  findings  of  this  search  are  given  in 
Appendix  A. 

2.  THE  AUTOMATED  CALIBRATION  CODE 

2.1  Purpose  and  Capabilities 

In  its  most  general  form,  the  bourtding  surface  model  requires  the  determination 
of  19  separate  constitutive  parameters,  including  2  initial  state  properties,  5  traditional 
material  constemts,  whose  values  may  be  directly  obtained  from  simple,  easy  to  perform 
laboratory  experiments,  and  12  model  constants,  which  must  be  indirectly  established 
through  a  trial  and  error  curve  fitting  process  using  the  results  of  undrained  triaxial 
tests.  A  general  summary  of  the  various  properties  is  presented  by  Herrmann  et  al 
(1980),  and  a  more  detailed  description  of  both  the  qualitative  and  quantitative  influence 
of  each  parameter  is  provided  by  OeNatale  (1982). 

This  breakdown  of  model  constants  is  common  to  most,  if  not  all,  of  the  soil 
model  formulations  introduced  in  recent  years.  Determination  of  the  directly 


measureable  or  "fixed”  parameters  is  straightforward  and  readily  accomplished. 
Determination  of  the  remaining  "free"  parameters,  however,  can  make  the  calibration 
procedure  prohibitively  difficult.  Rather  than  being  measured  directly  from  a  particular 
portion  of  a  specific  laboratory  test,  these  so-called  "free"  parameters  must  be 
determined  by  trial  and  error,  with  the  objective  being  to  obtain  the  best  overall  fit 
to  a  given  experimental  relation  or  set  of  observed  responses.  As  a  result,  the  overall 
accuracy  and  efficiency  of  the  calibration  process  can  be  strongly  dependent  on  both 
the  subjectivity  of  the  user  as  well  as  his  expertise  with  the  particular  material  model. 

In  formulations  such  as  the  bounding  surface  model,  which  employ  a  small  number 
of  material  parameters  whose  roles  in  the  constitutive  formulation  are  each  well 
defined,  the  manual  calibration  process  becomes  systematic  and  straightforward. 
However,  reliance  on  user  expertise  is  still  high,  since  all  manual  curve  fitting 
procedures,  by  their  very  nature,  require  both  judgement  (in  deciding  just  what 
constitutes  the  "best"  overall  fit)  and  familarity  (in  deciding  how  much  each  parameter's 
value  must  be  changed  to  improve  a  given  prediction). 

In  order  to  simplify  the  model  calibration  process,  a  computer  aided  calibration 
procedure  has  been  developed  and  tested.  Since  the  calibration  of  a  material  model 
involves  minimizing  the  error,  or  residual,  between  the  observed  and  predicted  soil 
response,  the  process  can  quite  naturally  be  viewed  as  an  optimization  problem.  Hence, 
the  computer  code  employs  a  quasi-Newton  optimization  strategy  to  locate  that  set 
of  parameter  values  which  minimizes  the  discrepancy  between  the  model  predictions 
and  the  experimental  observations  included  in  the  calibration  data  base.  The  code 
permits  any  number  of  tests,  relations  and/or  individual  observations  to  be  included  in 
the  calibration  data  base.  Different  weights  may  be  assigned  to  specific  components 
of  the  data  base  if  it  is  felt  that  certain  tests,  relations  or  observations  are  more 
reliable  or  representative  than  others,  or  if  it  is  necessary  to  have  the  final  model 
predictions  fit  some  data  more  closely  than  others.  Because  this  new  computer  aided 
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procedure  greatly  reduces  the  dependence  of  calibration  success  on  user  expertise,  it 
significantly  increases  the  accessibiiity  and  usefulness  of  sophisticated  material  models 
to  the  general  engineering  community.  Although  the  code  was  developed  specifically 
for  use  with  the  bounding  surface  plasticity  model,  it  can  be  readily  adapted  to  ot^er 
constitutive  formulations. 

In  addition  to  providing  a  means  for  automatically  calibrating  the  bourtding 
surface  model,  the  code  can  also  be  used  to  generate  a  set  of  model  predictions  for 
homogeneous  test  conditions.  The  code  may  thus  be  used  in  those  applications  where 
the  driving  program  EVAL  and  subroutine  CLAY  would  formerly  have  been  employed, 
Herrmann  et  al  (1980).  A  comprehensive  discussion  of  the  new  code,  including  a 
comparison  of  the  effectiveness  of  the  manual  and  automated  calibration  procedures, 
is  presented  by  OeNatale  (1982). 


2.2  The  Calibration  Data  Base 

The  ultimate  goal  of  the  calibration  process  is  to  identify  that  set  of  parameter 
values  which  enables  the  theoretical  model  to  most  closely  simulate  the  observed 
material  response.  This  goal  is  ordinarily  accomplished  by  fitting  the  model  to  a 
representative  set  of  experimentally  observed  stress>strain  relations  or  "calibration  data 
base".  Ideally,  this  calibration  data  base  should  be  complete  and  diverse  enough  that 
all  important  aspects  of  the  material's  response  are  included,  and  all  necessary 
constitutive  parameters  may  be  uniquely  established. 

In  its  most  general  form,  the  bounding  surface  formulation  is  a  fully  three- 
dimensional  stress-strain  model.  With  a  single  set  of  parameter  values  the  model  may 
be  applied  to  soil  stress  states  at  all  overconsolidation  ratios  (OCR's),  in  either 
compression  or  extension,  under  both  drained  and  undrained  conditions,  and  for  both 
monotonic  and  cyclic  loading.  Hence,  to  establish  the  optimal  values  of  the  necessary 
constitutive  parameters,  the  calibration  data  base  should  ideally  contain  observations 
from  the  following  seven  standard  laboratory  tests: 


1)  an  isotropic  (or  K^)  consolidation  or  drained  compression  test  with 
both  loading  and  unloading;  and, 

2-7)  undrained  triaxial  compression  and  extension  tests  on  specimens  in 
the  normally  (OCR  =  1),  Ughtly  (1  <OCR  <  2.5)  and  heavUy  (OCR  >  <f) 
overconsoiidated  ranges. 

The  results  of  the  consolidation  test  are  required  to  establish  the  slopes  of  the  isotropic 
consolidation  and  swell/recompression  curves  (in  e  -  In  p*  space)  X  and  ic.  These  two 
parameters  are  traditional  soil  properties  and  would  normally  be  assigned  values 
immediately,  prior  to  using  the  automated  calibration  procedure.  The  results  of  the 
six  undrained  triaxial  experiments  are  required  to  determine  the  12  model  constants 
cited  in  section  2.1,  and  would  thus  provide  the  data  needed  to  direct  the  automated 
calibration  procedure. 

The  triaxial  test  results  should  ideally  be  represented  in  terms  of  the  observed 
q  vs  p',  q  vs  and  u  vs  Cj  relations.  Naturally,  if  a  less  general  form  of  the  bounding 
surface  model  is  acceptable,  the  number  of  constitutive  parameters  involved  and  the 
number  of  laboratory  experiments  required  can  be  dirastically  reduced.  For  example, 
if  the  model  is  only  to  be  applied  to  normally  consolidated  soils  loaded  in  compression, 
the  number  of  required  constitutive  parameters  drops  from  19  to  7,  and  only  the 
isotropic  consolidation  and  a  single  triaxial  test  are  needed  for  model  calibration. 

Although  the  above  data  base  is  recommended,  the  bounding  surface  model  could 
also  be  calibrated  using  other  types  of  data.  For  example,  drained  rather  than  undrained 
tests  could  be  employed.  However,  untb’ained  tests  are  preferable,  sirxre  good  initial 
estimates  of  many  of  the  model  parameters  can  be  made  by  examining  the  expe-imentaliy 
observed  uxlrained  stress  paths. 

There  is  also  some  evidence  that  the  calibration  data  base  need  not  necessarily 
include  data  from  aU  three  overconsolidation  re^ons  (see  DeNataie  19S2).  That  is,  it 
may  be  sufficient  to  include  only  tests  from  Ite  normal  and  heavy  ranges,  or  perhaps 
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even  the  heavy  range  alone.  The  data  which  supports  this  poss  bility  is  not,  however, 
conclusive,  and  therefore  testing  at  all  three  regions  is  still  recommended. 

In  addition,  the  experimental  observations  need  not  necessarily  include  all  three 
relations  q  vs  p',  q  vs  Cj  and  u  vs  Cj.  Of  the  four  undrained  response  parameters  q, 
p',  u  and  only  three  are  independent.  In  practice,  p'  is  never  actually  measured, 
but  rather  is  computed  from  the  relation: 

3  (o,  -  u)  +  q 

p.  =  - 3 - 

3 

Hence,  any  two  of  the  three  relations  cited  above  will  completely  define  the  soil 
response.  The  use  of  q  vs  p'  or  q  vs  data  only  is  insufficient,  since  each  of  these 
relations  is  insensitive  to  certain  of  the  constitutive  parameters.  There  is  some 
evidence  that  the  use  of  u  vs  data  alone  may  be  sufficient  (see  DeNataie  1982). 
However,  the  use  of  ail  three  response  relations  appears  to  increase  the  rapidity  with 
which  the  optimization  algorithm  converges  to  the  minimum.  Presumably,  the  inclusion 
of  redundant  data  reinforces  the  correct  search  direction.  Since  the  cost  of  an 
automatic  calibration  run  is  only  mar^nally  affected  by  the  numbers  of  response 
.ations  included  in  the  calibration  data  base,  it  is  recommended  that  all  three  of 
the  relations  cited  above  be  used. 

Finally,  it  may  be  possible  to  use  other  than  triaxial  tests  to  acquire  the 
necessary  experimental  information.  Although  the  conventional  triaxial  apparatus  is 
the  most  common  and  versatile  laboratory  device,  the  simple  shear  apparatus  could 
also,  for  example,  be  used.  In  general,  the  soils  observed  stress-strain  characteristics 
will  be  to  some  extent  dependent  on  the  testing  device  which  is  used.  Thus,  in  practial 
problems,  the  laboratory  device  used  to  generate  the  calibration  data  base  should 
simulate,  as  closely  as  possible,  the  loading  conditions  for  which  the  bounding  surface 
predictions  will  eventually  be  generated. 
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2.3  Practical  Considerations 

The  goal  of  the  calibration  procedure  is  to  identify  that  set  of  parameter  values 
which  minimizes  the  error  between  the  model  predictions  and  the  experimental 
observations  included  in  the  calibration  data  base.  The  automated  optimization  code 
thus  seeks  to  locate  the  objective  function’s  global  minimum.  Unfortunately,  there  is 
no  guarantee  that  the  algorithm  will  always  succeed.  The  quasi>Newton  strategy 
employed  by  the  model  calibration  code,  like  most,  if  not  all,  practical  optimization 
algorithms,  is  designed  only  to  locate  local  minima  in  the  vicinity  of  the  initial 
estimates.  Hence,  the  probability  that  the  true  global  minimum  will  be  found  is 
directly  related  to  the  degree  of  unimodality  exhibited  by  the  objective  function  and 
the  accuracy  of  the  initial  guess. 

Preliminary  research  by  DeNatale  (19S2)  indicates  that  the  use  of  the  "absolute 
-  Euclidean"  measure  of  error  leads  to  a  more  unimodal,  and  thus  desireable,  objective 
fixiction.  A  procedure  for  acquiring  improved  starting  estimates  has  also  been  developed 
by  DeNatale  (19S2).  Through  actual  testing  with  a  number  of  different  soils,  this 
strategy  has  been  found  to  produce  starting  estimates  which  enable  the  automated 
calibration  code  to  consistently  locate  the  global  minimum.  In  practice,  however,  the 
only  way  to  ensure  that  the  global  minimum  has  been  found  is  to  conduct  the  search 
from  a  variety  of  different  starting  points.  The  solution  which  yields  the  lowest  value 
of  the  objective  function  may  then  be  regarded  as  the  global  minimum. 

A  second  practical  consideration  concerns  the  quality  of  the  calibration  data 
base.  The  user  should  ensure  that  the  experimental  observations  included  in  the 
calibration  data  base  are  diverse  enough  to  permit  the  optimal  values  of  the  required 
ixiknown  model  parameters  to  be  iniquely  defined.  For  example,  if  the  code  is  used 
to  identify  those  model  parameters  associated  with  heavily  overconsoii dated  material 
response,  the  calibration  data  base  must  include  observations  made  on  heavily 
overconsolidated  specimens.  If  the  necessary  experimental  data  is  not  included,  the 
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program  will  continue  to  execute,  but  the  final  computed  "optimal"  values  of  the 
undefined  parameters  will  be  very  close  to  the  initial  estimates.  The  major  consequence 
of  an  inadequate  or  incomplete  calibration  data  base  is  related  to  the  cost  of  the 
analysis.  Certain  computational  costs  increase  in  proportion  to  n  (n  =  number  of 
parameters  to  be  determined),  and  a  single  gradient  evaluation  requires  either  n  or  2n 
additional  objective  furxrtion  evaluatiorts  (depending  on  whether  forward  or  central 
differencing  formulae  are  used).  Thus,  to  minimize  the  cost  of  a  model  calibration 
run,  the  user  should  seek  to  identify  only  those  parameters  whose  values  can  be  defined, 
given  the  particular  data  base.  A  comprehensive  discussion  of  the  influence  of  each 
of  the  19  model  parameters  is  provided  by  DeNatale  (1982),  arwi  may  be  referred  to 
if  uncertainly  exists. 

2.#  Example 

To  verify  the  viability  of  the  new  computer  aided  calibration  procedure,  the 
method  was  applied  to  a  number  of  representative  data  bases,  both  artificial  and  real. 
The  outcome  of  these  studies  is  discussed  by  DeNatale  (1982).  Among  the  real  data 
bases  to  which  the  automated  process  was  applied  are  the  experimental  results  on 
Kaolin  reported  by  3afroudi  (1982). 

The  Bounding  Surface  model  was  calibrated  on  the  basis  of  conventional  undrained 
triaxial  compression  and  extension  tests  on  samples  at  overconsolidation  ratios  of 
OCR  =  1,  2  and  6.  With  the  necessary  constitutive  parameters  thus  having  been  fixed 
at  their  optimal  values,  predictions  were  then  generated  for  a  variety  of  additional 
undrained,  drained,  hollow  cylinder  and  cyclic  tests. 

The  results  obtained  with  the  automated  procedure  are  compared  to  those  of 
the  manual  solution,  as  reported  by  3afroudi  (1982)  and  Herrmann  et  al  (1982),  in 
Table  1  and  Figures  1  through  6.  As  may  be  observed  in  Table  1,  the  optimal  values 
of  the  model  parameters  as  established  through  the  automated  and  manual  procedure 
are,  as  a  group,  distinctly  different.  For  the  given  set  of  options  specified  in  Table 
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1,  the  automated  solution  was  better  than  the  hand  solution,  in  the  sense  of  having  a 
lower  objective  function  value.  The  relative  merits  of  the  two  solutions  may  perhaps 
most  clearly  be  seen  by  comparing  the  associated  calibration  curves  shown  in  Figures 
1  through  6.  On  the  whole,  the  automated  curves  appear  to  more  closely  match  the 
experimental  observations.  A  point  by  point  comparison  of  the  two  solutions  is 
essentially  impossible,  since  in  performing  a  manual  calibration,  different  sifbjective 
weights  are  implicity  assigned  to  the  various  components  of  the  data  base  whidi  cannot 
be  precisely  identified. 

2.5  Cost 

The  automated  calibration  code  has  been  written  in  FORTRAN  and  implemented 
or  both  an  LSI-11/23  minicomputer  as  well  as  a  VAX-11/780  super-minicomputer.  The 
cost  of  a  given  analysis  is  controlled  primarily  by  the  number  of  distinct  experimental 
tests  included  in  the  calibration  data  base  and  the  number  of  constitutive  parameters 
whose  optimal  values  are  being  sought.  A  typical  computer  calibration,  such  as  the 
11-dimensional  problem  reported  herein  for  the  data  of  Jafroudi  (1982),  requires  from 
400-600  objective  function  evaluations,  or  about  60-90  minutes  of  VAX  CPU  time,  at 
a  cost  of  approximately  $50.00-$75.00.  This  cost  is  relatively  low  when  compared  to 
the  cost  of  the  experimental  program  needed  to  establish  the  calibration  data  base, 
and  in  li^t  of  the  resulting  economy  of  finite  element  analyses  based  on  the  model. 
Calibration,  of  course,  need  only  be  done  once  for  a  given  soil,  regardless  of  the 
rxjmber  and  variety  of  svi^sequent  finite  eiement  analyses  which  utilize  the  model. 

3.  NUMERICAL  IMPLEMENTATION 

3.1  bicrementalization  Of  Rate  Equations 

Using  tensor  rotation,  the  basic  rate  equations  for  the  bounding  surface  plasticity 
model  can  be  written  in  the  form  (repeated  indices  are  summed  from  I  3  and  free 
indices  take  on  values  of  1,  2  or  3.) 
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(1) 


The  effective  stress  components  are  denoted  by  o..  and  the  strain  components  by  e-|. 
In  general,  D...  ,  is  a  function  of  both  o..  and  C;;  and  one  or  more  internal  variables 

IJKI.  1)  1) 

q..  The  specific  form  of  for  the  bounding  surface  model  for  cohesive  soils  is 

given  by  Herrmann  et  ai  (1980). 

For  numerical  analysis  purposes,  it  is  more  convenient  to  express  the  relationship 


in  matrix  form;  i.e.. 


{o}  =  IDI  {e} 


where  ([  ]  '  is  the  matrix  transpose) 

(<4'^  =  (<J,,  Oy.  V  T„,  T^) 

{clT  =  (c^,  ty.  e^,  1^,  1f„,  (3) 

The  tensor  components  of  shear  strain  e..  are  one>half  of  the  engineering 
components  y„.  The  symmetric  6x6  iDl  matrix  is  expressed  in  terms  of  the  components 
of  the  3x3x3x3  tensor  by  using  the  following  six  sets  of  corresponding  indices 

(1;1,1)  (2;2,2),  (3;3,3),  (4;1,2),  (5;1,3)  and  (6;2,3)  where  the  first  number  is  the  row 
(column)  number  in  the  matrix  and  the  second  and  third  numbers  are  the  first  Oast) 
two  indices  for  the  tensor. 

To  be  able  to  use  eq.  (2)  in  an  incremental  solution  procedure,  it  must  be 
expressed  in  an  incremental  form.  Consider  the  n'*'  step  of  an  incremental  analysis; 
i.e.,  the  solution  has  been  found  at  n-1,  and  it  is  desired  to  calculate  the  incremental 
change  that  will  give  the  solution  at  n.  Because  of  the  nonlinear  behavior,  iteration 
is  required  to  establish  the  ir^remental  change.  In  the  k-1*^  iteration  of  this  process, 
the  estimates  of  the  stress  and  strain  states  at  n  are  given  by 

('•n,k.l  =  ‘  <'*>n.k-l 
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Even  though  rate  independent  behavior  is  being  considered,  it  is  convenient  to 
think  in  terms  of  the  time  histories  of  the  quantities  involved.  Integrating  eq.  (2) 
from  time  t^j  to  t^^,  and  using  the  trapezoidal  formula  to  approximate  the  right  hand 
side,  gives 


(6) 

®n.k-l  *  ‘“J.k-l 

(7) 

When  eq.  (6)  is  used  in  a  finite  element  analysis  it  is  tacitly  assumed  that  all 
the  stress  and  strain  components  at  a  particular  point  in  the  body,  and  from  point  to 
point,  change  proportionally  from  their  values  at  n-1  to  their  values  at  n.  Thus,  in 
order  that  the  true  solution  history  be  accurately  modeled,  as  required  for  an  inelastic 
material,  the  solution  step  size  must  be  limited. 

Eq.  (6)  is  the  desired  incremental  stress-strain  equation  for  iteration  k  of 
increment  n.  Because  [D]^  is  a  function  of  the  stress  and  strain  states  at  n  (due  to 
the  dependence  of  on  o.j  at  Cjj),  it  is  necessary  to  base  its  value  on  the  estimates 
of  the  previous  iteration  (eqs.  (4)  and  (5)).  The  resulting  value,  denoted  by  [D]  .  ,, 
is  used  in  eq.  (7).  The  fact  that  is  also  a  function  of  the  stress  and  strain 

states  at  n  is  less  obvious.  This  dependence  arises  because  of  the  difference  in  material 
response  for  loading  and  unloading  conditions  and  the  fact  that  whether  loading  or 
unloading  occurs  during  a  given  increment  is  influence  by  the  values  of  and  {c}^. 
That  is,  lDlp_2  is  the  tangent  stiffness  at  the  beginning  of  the  increment,  and  this 
stiffness  differs  for  loading  and  unloading  conditions,  as  determined  by  the  values  of 
and  {Ae}^.  The  appropriate  value  of  is  written  as  j. 

At  the  beginning  of  the  iteration  process,  initial  estimates  are  required.  For 
the  first  iteration  ot  the  first  irx:rement,  they  are  usually  taken  to  be  zero.  For  the 
initial  iterations  of  succeeding  increments,  they  also  can  be  started  at  zero;  however. 
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it  may  be  desirable  to  make  use  of  information  from  the  previous  increment  to  obtain 

better  starting  values.  The  simplest  procedure  is  to  use  as  the  initial  estimate  tHe 

final  values  fourKl  in  the  previous  irKrrement.  This  practice  is  based  on  the  assumption 

of  relatively  isiiform  behavior  from  irxrrement  to  increment.  Difficulties  can  arise 

when  the  histories  of  the  applied  external  agents  acting  on  the  structure  cause  a  switch 

from  loading  to  unloading  in  an  urutable  material  resporrse  regime,  and  these  external 

agents  are  loads  not  displacements.  For  example,  consider  the  one-dimensional  response 

shown  in  Figure  7.  Consider  the  case  when  the  state  of  the  soil  is  at  point  "A"  at 

the  end  of  increment  n-1.  If  during  increment  n,  is  specified,  two  final  states 

B  and  B*  are  possible.  One  corresporxis  to  Ac  (negative)  and  the  other  to  Ae*  (positive). 

Without  any  additional  information,  no  choice  can  be  made  between  B  smd  B*.  (It  is 

easily  seen  that  for  specified  stress  increments  in  the  stable  region  of  behavior  and 

for  specified  strain  inctements  anywhere,  no  such  problem  exists.)  The  suggested 

solution  to  this  impasse  is  to  assume  that  the  user  would  not  attempt  a  stress  controlled 

specification  for  "loading"  conditions  (path  A-B')  in  an  unstable  region  and,  hence,  if 

the  stress  increment  is  specified,  unloading  is  the  proper  behavior  (path  A-B)^.  For 

stress  controlled  conditions,  the  selection  of  the  unloading  path  can  be  assured  if  the 

starting  estimate  of  strain  is  of  opposite  sign  to  that  calculated  in  the  previous 

increment.  Thus,  the  following  strategy  is  recommended.  When  considering  a  series 

of  increments  for  which  the  rates  of  the  externally  applied  loads  and  displacements 

do  not  change  sign  ,  { Ae}  and  {Ac}  are  used  as  starting  estimates  for  increment 

n-1  n-1 

^  This  argument  requires  that  the  arrival  at  A  must  have  been  preceeded  by  strain 

controlled  steps. 

It  is  assumed  that  this  condition  is  sufficient  to  prevent  a  general  switch  from 


loading  to  unloading  within  the  soil  mass. 


n.  However,  as  one  such  solution  history  segment  is  ended  and  a  new  one  begins,  the 
conditions  necessary  for  the  non-uniqueness  problem  may  occur.  Hence,  for  the  first 
increment  of  each  such  series,  it  is  suggested  that  the  starting  strain  estimate  be 
taken  as  some  small  negative  multiple  (e.g.  -.01)  of  the  value  found  in  the  previous 
increment  (the  stress  increment  would  be  used  tnchanged).  The  reduction  in  absolute 
magnitude  is  in  deference  to  the  greater  stiffness  encountered  in  unloading.  Such  an 
initial  estimate  will  force  the  solution  to  select  path  A  B  if  the  necessary  conditions 
exist  for  the  above  mentioned  non-tsiiqueness  to  occur.  If  non-isiiqueness  is  not  a 
problem,  the  only  effect  of  this  procedure  is  to  slightly  slow  the  convergence  process. 

It  is  important  to  note  that,  in  general,  the  estimates  of  the  stress  and  strain 
increments  used  in  the  calculation  of  ^  do  not  in  fact  satisfy  eq.  (6).  The 

consequences  of  this  inconsistency  will  be  discussed  later. 

A  FORTRAN  subroutine  CLAY  for  the  calculation  of  the  matrix  has 

been  written  and  is  documented  by  Herrmann,  et  al  (19S1). 

3.2  Calculation  Of  Pore  Water  Pressure 

The  bounding  surface  plasticity  theory  is  expressed  in  terms  of  effective  stress, 
whereas  most  soil  related  problems  involve  the  application  and  calculation  of  total 
stress.  The  total  and  effective  stresses  differ  by  the  pore  water  pressure  u.  There 
are  three  posabilities  concerning  the  development  of  pore  water  pressure  in  soil:  ideal 
drained  conditions  (where  the  pore  water  pressure  is  identically  zero),  ideal  undrained 
conditions  (where  the  soil  is  completely  saturated,  and  no  flow  of  water  occurs),  and 
the  more  realistic  situation  where  there  is  a  global  flow  of  water  and/or  the  filling 
of  voids.  In  many  analyses  ideal  drained  or  undrained  conditions  are  assumed,  even 
though  they  may  only  be  approximately  true. 

The  total  stress  increment  is  the  sum  of  the  effective  stress  increment  and 
the  pore  water  pressure  increment: 
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For  drained  conditions  u=0  and  and  eq.  (6)  is  the  desired  relationship 

between  the  total  stress  and  strain  increments. 

For  undrained  conditions  there  are  several  possible  ways  of  proceeding.  The 
traditional  approach  has  been  to  neglect  the  (slight)  compressibility  of  the  water  and 
the  soil  particles,  and  thus  assume  incompressible  material  behavior.  However,  the 
finite  element  analysis  of  incompressible  materials  requires  a  special  formulation 
(Herrmann  1965;  Zienkiewicz  1977). 

In  order  to  avoid  having  to  deal  with  separate  formulations  for  drained  and 
undrained  conditions,  it  is  convenient  to  express  them  in  a  common  form.  This  can 
be  accomplished  if  the  slight  compressibility  of  the  soil  particles  and  the  pore  water 
is  recognized,  Sangrey,  et  al  (1969).  (An  alternative  interpretation  is  to  consider  the 
undrained  soil  as  incompressible,  and  to  approximately  specify  the  condition  by  means 
of  a  “penalty  function^',  Zienkiewicz,  et  al  1981,  where  the  associated  "penalty  number" 
corresponds  to  the  bulk  modulus  D.  'Rius,  the  pore  water  pressure  u  is  written  in 
terms  of  the  combined  bulk  modulus  F  of  the  soil  particles  and  the  pore  water,  and 
the  resulting  (very  small)  volume  change  Cj^j^  i.e.,  u=r  c^.  As  the  soil  becomes 
incompressible.  Drained  conditions  are  obtained  when  r=0.  For  indrained  conditions 
the  value  of  F  is  very  large  compared  to  the  terms  in  [D]  .  Thus,  the  soil  behaves 
as  a  "nearly  incompressible  soli<f'  (Herrmann  1965)  and  care  must  be  exercised  to  avoid 
numerical  round-off  and  excessive  constraint  problems.  Two  approaches  are  commonly 
used  to  achieve  this  goal.  One  method  is  to  use  the  special  formulation  given  by 
Herrmann  (1965)  for  incompressible  and  nearly  incompressible  solids,  while  the  other 
is  to  use  "reduced"  or  "selective-reduced"  integration  (Zienkiewicz  1977)  for  the  element 
stiffness  matrix;  the  importance  of  selecting  a  proper  element  type,  in  order  to  achieve 
acceptable  accuracy,  is  discussed  by  Nagtegaal,  et  al  (1974)  and  Zienkiewicz,  et  al 
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(1981).  In  the  latter  case,  the  above  expression  is  used  to  eliminate  u  from  eq.  (8); 
i.e.. 


•t  •  ».  •  » 

o. .  =  o. .  +  r  e, .  o. . 
»)  kk  ij 


(9) 


Integration  over  increment  n  ^ves: 


AcjT.  =  to.-  ♦  r  Acaj  • 

"n  "n  “ 

Eliminating  Aa..  using  eq.  (6)  and  retinming  to  matrix  notation  yields: 


(10) 


(11) 


where 


®’ln,k.l  '  ®n,k-l  *  '«  <■« 

All  components  of  [dl  are  zero,  except  i 

Returning  to  the  more  realistic  situation  where  water  movement  takes  place, 
two  cases  can  be  distinguished.  The  first  occurs  when  there  b  no  global  movement 
of  water  (either  the  time  scale  is  too  short  for  significant  flow  to  occur,  or  the  soil 
is  stressed  homogeneously,  thus  producing  no  pressure  gradients)  in  a  partially  saturated 
soil.  This  condition  can  be  modeled  by  using  a  variable  F  which  is  a  function  of  the 
current  saturation  state. 

When  there  is  actual  global  flow  of  water,  it  is  necessary  to  perform  a  coupled 
flow-stress  analysis  (Sandhu  and  Wilson  1969).  The  details  of  sudi  an  analysis  are 
beyond  the  scope  of  this  study.  It  should  be  noted  that  the  bounding  surface  plasticity 
model  is  valid  for  such  situations;  however,  the  pore  water  pressure  can  not  simply 
be  calculated  from  tfie  expression  Fcj^j^.  Eq.  (6)  is,  however,  still  valid  for  relating 
the  increments  of  effective  stress  and  strain. 
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3.3  General  Discuuion  Of  Finite  Element  Implementation  Of  Plasticity  Models 

The  desired  characteristics  of  any  numerical  scheme  are  ease  of  implementation, 
computational  efficiency  and  robustness  of  behavior.  It  is  on  the  bases  of  these  three 
characteristics  that  the  competing  methods  will  be  judged.  There  is  of  course  a 
considerable  body  of  literature  available  on  the  analysis  of  elasto-plastic  bodies.  Only 
those  items  which  directly  relate  to  this  work  will  be  discussed  here;  for  a  more 
general  discussion  the  reader  is  referred  to  Owen  and  Hinton  (1980)  and  Nagtegaal 
et  al  (1974). 

In  general,  the  response  of  an  elasto-plastic  body  is  highly  nonlinear  and  path 
dependent.  Thus,  a  general  numerical  analysis  procedure  for  elasto-plasticity  problems 
requires  an  incremental  solution.  Unless  the  increments  are  made  excessively  small, 
iteration  must  be  conducted  in  each  increment  to  account  for  the  nonlinear  behavior. 
The  two  most  commonly  employed  classes  of  iterative  methods  are  successive 
approximation  (substitution)  and  Newton  like  procedures.  The  method  of  successive 
approximation  can  be  cast  in  a  variety  of  forms  and  thus  is  not  a  utique  operation 
(see  Isaacson  and  Keller  1966).  The  alternative  forms  range  from  extremely  simple, 
but  very  slowly  convergent,  procedures,  to  more  complicated  methods.  The  Newton- 
Raphson  method  and  numerous  approximations  to  it  constitute  the  class  of  "Newton 
like"  methods.  The  many  available  solution  methods  obviously  could  not  all  be  evaluated 
in  this  study.  The  selection  of  the  methods  that  are  compared  in  this  study  is  discussed 
in  the  next  section. 

3.%  Selection  Of  Methods  For  Comparison 

While  the  method  of  successive  approximation  is  extremely  easy  to  implement, 
it  often  suffers  from  poor  convergence  characteristics.  For  nonlinear  elasticity 
problems,  where  the  entire  solution  history  can  be  accounted  for  in  one  step,  there 
appears  to  be  little  question,  unless  the  non  linearities  are  very  weak,  that  a  Newton 
like  method  is  to  be  preferred.  However,  for  inelasticity  problems  where  relatively 


small  steps  are  required  to  account  for  the  path  dependence  of  the  solution,  the  choice 
is  not  so  clear.  The  question  is  further  clouded  by  the  fact  that,  for  a  given  solution 
step,  and  at  a  given  point  in  an  elasto-plastic  body,  the  stiffness  may  be  disoontinuous 
due  to  the  onset  of  yielding  or  to  the  progression  from  (local)  loading  to  unloading 
conditions.  The  latter  case  is  of  concern  in  highly  statically  indeterminate  situations 
where,  at  a  given  point  in  space  and  time,  it  is  not  known  a  priori  whether  or  not 
the  material  will  experience  loading  or  tmloading.  This  stiffness  discontinuity  may 
make  convergence  of  Newton  like  method  slower  than  for  problems  where  the  Sacobian 
is  strictly  continuous. 

In  order  to  keep  the  scope  of  dw  study  within  practical  bounds,  certain 
acceptability  criteria  are  stated  aixi  used  to  limit  the  number  of  methods  to  be 
compared.  Two  of  the  key  reasons  for  the  computational  efficiency  of  the  finite 
element  method,  as  applied  to  structural  problems,  are  the  symmetry  and  baiMied  nature 
of  the  iumultaneous  equations.  Hence,  for  this  study,  it  is  required  that  the  solution 
methods  preserve  these  characteristics. 

It  is  on  the  basis  of  maintaining  symmetry  that  Owen  and  Hinton  (1980)  rule 
out  the  general  form  of  the  Newton-Raphson  method,  arK)  instead  advocate  the  use  of 
an  approximate  form;  i.e.  the  "tangent  stiffness”  method. 

Currently  there  is  considerable  interest  in  quasi-Newton  methods  for  nonlinear 
structural  problems  (see  Geradin  et  al  1981).  The  central  goal  of  qua^-Newton  methods 
is  to  avoid  calculating  the  Jacobian  every  iteration.  Instead,  simple  tpdating  formulas 
are  used  to  approximate  the  Jacobian  (or  its  inverse)  in  terms  of  the  previous 
approximation  and  simple  vector  quantities  (the  previous  solution  and  current  residual). 
The  ipdates  can  be  made  directly  to  the  Jacobian  or  to  its  inverse;  in  the  first  case, 
a  set  of  simultaneous  equations  must  be  solved  each  iteration,  while  in  the  second, 
only  a  matrix  multiplication  is  required,  bt  the  former  case,  because  of  the  continued 
need  of  solving  a  set  of  simultaneous  equations  at  each  iteration,  arxl  because  the  cost 


of  using  the  ipdating  formula  is  as  great  as  the  cost  of  calculating  the  approximate 
3acobian  for  the  other  methods  considered  in  this  study,  there  would  seem  to  be  little 
or  no  advantage  offered  by  the  method  (although  in  optimization  and  nonlinear  elasticity 
problems,  other  important  advantages  make  it  a  viable  method).  Thus,  only  those 
quasi-Newton  methods  which  ig>date  the  inverse  would  appear  to  be  of  interest  in 
elasto-plasticity  problems.  Inverse  igadating  methods,  however,  result  in  dealing  with 
a  full  matrix,  thus  destroying  the  advwitages  of  the  banded  characteristic  of  finite 
element  equations.  The  ideal  situation  would  be  an  updating  formula  which  can  be 
applied  directly  to  the  reduced  (tpper  triangular  form)  Jacobian  and  which  would  not 
destroy  its  banded  nature  nor  its  original  symmetry  (a  characteristic  required  for  the 
efficient  reduction,  each  iteration,  of  the  load  vector).  The  available  quasi-Newton 
schemes  of  this  type  do  not  satisfy  the  ease  of  implementation  requirement  (see  the 
discussions  by  Geradin  et  al  19S0,  1981,  Mathies  and  Strang  1979  and  Schubert  (1970). 
Thus,  a  simpler  form  is  considered  herein.  Finally  it  should  be  noted  that  because  of 
the  lack  of  a  natural  objective  fisKtion,  the  line  search  criterion  of  the  quasi-Newton 
methods,  as  applied  to  optimization  (see  Fletcher  1980)  and  nonlinear  elasticity,  is  lost 
for  inelasticity  problems.  While  alternative  criteria  have  been  proposed  for  elasto- 
plasticity  problems,  they  appear  to  lack  the  simplicity  and  robustness  of  the  minimization 
of  an  objective  function. 

The  chief  appeal  of  inverse  quasi-Newton  methods  is  the  elimination  of  the 
necessity  of  reducing  the  Jacobian  every  iteration.  This  same  objective  can  be  achieved 
by  using  a  modified  iteration  method  which  only  occasionally  updates  the  Jacobian. 
At  other  times  the  effects  of  the  changes  in  the  Jacobian  are  estimated  and  transferred 
to  the  right-hand  side  of  the  equations.  This  procedure  is  highly  recommended  by 
Owen  and  Hinton  (1980)  and  is  included  in  this  evaluation. 

Successive  approximation  methods  are  popular  because  of  their  ^mplicity  of 
implementation.  Thus  one  form  of  successive  approximation  is  studied,  h  order  to 
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demonstrate  its  relation  to  the  Newton>Raphaon  method  it  is  derived  as  an  alternative 
to  the  tangent  stiffness  method. 

bi  order  to  improve  the  rate  (and  also  to  enlarge  the  domain)  of  convergence 
of  the  iteration  process,  acceleration  schemes  are  often  used.  These  schemes  usually 
employ  some  type  of  extrapolation  in  order  to  obtain  a  better  solution  estimate  than 
given  directly  by  the  iteration  process  (see  Isaacson  and  Keller  1966).  The  extrapolation 
can  usually  be  epxressed  in  terms  of  an  acceleration  (iteration,  relaxation)  factor.  The 
simplest  method  employs  a  constant  factor  selected  by  the  user  on  the  basis  of  past 
experience.  Because,  for  inelasticity  analyses,  the  optimum  factor  can  vary  widely 
from  problem  to  problem,  and  even  from  increment  to  increment,  some  type  of  strategy 
for  automatically  selecting  it  is  desirable.  For  one  degree  of  freedom  problems,  the 
Aitken's  V  method  (see  Isaacson  and  Keller  1966)  is  simple  to  apply  and  has  proven 
to  be  quite  effective.  Two  methods  for  adapting  it  to  multi-degree  of  freedom  problems 
are  considered  herein. 

The  final  question  that  is  addressed  is  how  best  to  handle  the  inconsistency 
(previously  mentioned)  between  the  estimates  {AcJ  .  ,  and  {Ae}  .  ,  used  to  calculate 
and  eq.  (6);  three  schemes  are  compared. 

3.5  Theory 

3.5.1  Solution  Methodb 

In  the  following  discussion  it  will  be  assumed  that  the  reader  is  familiar  with 
the  standard  steps  involved  in  formulating  a  finite  element  analysis. 

For  a  given  solution  increment  eq.  (6)  is  used  in  the  formulation  of  element 
stiffness  matrices,  which  are  in  turn  combined  to  form  the  system  stiffness  matrix 
The  incremental  load  vector  is  denoted  by  Equilibrium  leads  to  the 

following  system  of  simultaneous  equations  for  the  incremental  displacements 

dCln^^ln  =  ^^Jn 
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Because  [iT]^  is  a  function  of  { and  { the  above  equations  are  nonlinear 
and  require  iteration  to  solve.  The  Newton-Raphson  method  gives  for  iteration  k: 


-  '’’nik-l  <‘*)n,k.l 


Where  the  residual  vector  is: 


<''*>n,k.l  =  I^'n.k-I  I'^ln.k-l  ’  <*’> 

The  components  of  the  3acobian  j  are  found  by  taking  derivatives  of  eq.  (15) 

with  respect  to  the  components  of  in  index  notation: 


-  F,)  r 

'‘i„.k-i '  — ®r~  n.k.i  “  hi  *  ® 


n,k-l 


'^n.k-l  =  »^n,k-l  •  i^iVk-l 

Owen  and  Hinton  (19S0)  state  that  in  general  {iD'  is  not  symmetric.  In  addition,  it 
is  relatively  difficult  and  e}q>ensive  to  compute.  Based  on  the  classical  graphical 
interpretation  of  the  Newton-Raphson  method,  Owen  and  Hinton  (19S0)  suggest  that, 
instead  of  using  eq.  (17),  the  Jacobian  be  approximated  by  the  tangent  stiffness  matrix 
at  "n”.  This  requires  using  only  the  j  matrices  (see  paragraph  following  eq. 

(7))  for  the  formation  of  the  system  tangent  stiffness  matrix,  call  the  result 
The  matrix  is  still  needed  for  the  calculation  of  the  residual  vector  eq.  (15). 

Thus  the  "tangent  stiffness  method"  consists  of  iteration  using  sqs.  (14)  and  (15)  with 


'^n.k.I  * 
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As  an  alternative,  eq.  (17)  can  be 


approximated  by  neglecting  the  last  term 


Inspecting  eq.  (21)  it  is  seen  that  this  second  approximation  to  the  Newton-Raphson 
method  is  one  of  the  possible  forms  of  the  classical  method  of  successive  approximations, 
and  shall  be  referred  to  by  that  name  in  the  remainder  of  the  report.  It  can  be 
applied  by  either  using  eq.  (21)  directly  or  using  eqs.  (14),  (15),  and  (19).  The  use  of 
equations  (14),  (15)  and  (19)  would  be  expected  to  be  somewhat  less  susceptible  to 
round-off  error;  however,  no  significant  differences  were  detected  in  the  examples 
analyzed  in  this  study. 

The  calculation  of  the  residial  vector,  eq.  (15),  requires  some  special  attention. 
One  option  is,  at  the  element  level,  to  use  j  to  calculate  an  clement  stiffness 

matrix,  to  then  multiply  this  matrix  by  the  k-1  estimate  of  the  displacements  for  the 
nodes  defining  the  eiement,  and  to  add  it  to  the  negative  of  the  element  load  matrix 
(i.e.,  a  direct  evaluation  of  eq.  (15)  at  the  element  level).  The  resulting  element 


^  When  convergence  occurs,  eq.  (13)  is  exactly  satisfied.  Thus,  neither  this 
approximation,  nor  the  one  leading  to  the  tangent  stiffness  method,  has  any  effect 
on  the  final  accuracy  of  the  solution. 
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residual  vectors  are  then  combined  in  the  usual  way  to  give  the  system  residual  vector. 
WhUe  this  operation  requires  little  additional  computational  effort  for  the  method  of 
successive  approximations,  it  does  for  the  tangent  stiffness  method.  The  cause  of  this 
additional  work  is  the  need  to  calculate  two  element  stiffness  matrices^  based  upon 
the  two  quantities  (D]^  and  l^_j.  The  first  element  stiffness  matrix  is  needed 
in  the  calculation  of  the  system  tangent  stiffness  matrix,  and  the  second  in  the 
calculation  of  the  residual  vector.  An  alternative  procedure  for  the  calculation  of  the 
residual  vector  avoids  this  additional  effort;  however,  it  places  certain  restrictions  on 
the  order  of  the  integration  used  in  establishing  the  element  matrices.  In  this  second 
approach,  an  initial  stress  vector  is  calculated: 

('“’oln.k-l  =  ®n,k.l  (‘^'n.k-l  <“> 

The  following  incremental  stress>strain  equation  is  then  used  in  the  calculation  of  the 
element  stiffness  and  load  matrices: 

f»‘*n,k  =  '“n.k-l  (‘'»n,k  ‘  (Hln.k-I  <”> 

Assuming  that  all  numerical  integrations  are  done  with  the  same  accuracy,  it 
is  easy  to  show  that  the  use  of  eq.  (23)  yields  element  matrices  that,  when  combined 
at  the  system  level,  ^ve  the  desired  tangent  stiffness  matrix  and  the  residual  force 
vector.  In  this  operation  care  must  be  taken  to  ensure  that  all  numerical  integrations 
are  of  the  same  order.  For  exampie,  for  a  four  node  element,  if  four  point  integration 
is  used  in  calculating  the  element  matrices  but  the  incremental  properties  (eq.  (23)) 
are  calculated  only  at  the  element  center  (i.e.,  assumed  constant  over  the  element), 

^  Geometric  quantities  such  as  the  shape  function  derivates  need  only  be  calculated 
once. 
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an  inconsistency^  exists  and  convergence  is  dirastically  affected.  Thus  when  using  this 
scheme  the  stresses,  strains  and  properties  must  be  evaluated  at  the  integration  points. 
However,  for  a  four  node  element  it  has  been  shown  (see  Herrmann  1972)  tint  the 
stress  calculations  are  most  accurate  at  the  element  center,  not  at  the  integration 
points.  For  this  reason,  and  because  the  four  node  element  is  poorly  behaved  for 
undrained  conditions,  the  use  of  an  eight  or  nine  rK>de  element  with  stresses  and 
properties  calculated  at  each  of  the  four  quadrature  points  is  recommended  for  future 
work. 

5.3.2.  Reduction  Of  The  Number  Of  Equation  Triangularizations 

The  use  of  the  "modifierf*  Newton  method  (referred  to  as  the  "initial  stiffness" 
method  by  Owen  and  Hinton  1980)  is  the  classical  means  for  reducing  the  number  of 
triangularizations  of  the  left>hand  ade  of  the  umultaneous  equations.  In  this  procedure, 
the  stiffness  matrix  is  updated  only  occasionally.  Because  of  the  drastic  difference 
in  stiffness,  in  elasto-plasticity  problems  (which  is  encountered  in  a  progression  from 
elastic  behavior  to  yielding  and  from  loading  to  inloading)  it  is  desirable  to  update  at 
least  once  each  increment.  For  this  study,  the  stiffness  matrix  b  updated  in  the 
second  iteration  of  each  increment.  The  second  iteration  was  chosen  because  it  often 
takes  at  least  one  iteration  to  establish  the  loading-unloading  characteristics  for  the 
increment.  In  addition,  it  is  igKiated  in  the  first  iteration  of  the  first  increment  of 
each  new  loading  segment  (a  loading  se^nent  generally  consists  of  many  increments). 
This  update  is  done  because  there  is  often  a  switch  from  loading  to  unloading  at  the 
beginning  of  a  new  loading  segment.  Finally,  the  stiffness  matrix  is  updated  every 
IRPET  iteration,  where  IRPET  is  specified  by  the  user.  The  use  of  the  "modified 
Newton"  procedure  in  conjunction  with  the  method  of  successive  approximations, 

^  That  is,  the  residial  node  point  forces  are  not  accurately  made  to  be  zero. 
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implemented  by  means  of  eq.  (21)  instead  of  eq.  (14),  requires  that  the  difference 
between  the  present  stiffness  matrix  and  the  last  stiffness  matrix  to  be  triangularized, 
be  accounted  for.  With  this  end  in  mind,  eq.  (21)  is  written  in  the  form: 


-  Ml)  ( 

Now  if  the  quantity  { Au)  i.  Of*  the  ri^t-hand  side  is  estimated  by  means  of  { Au}  , 

n|K  n^r 

(as  convergence  occurs  no  approximation  is  introduced): 


n,k 


(25) 


where 


(“’tS.k.l  =  (""In  -  {'“n.k-l  -  <‘“>„,k.| 


(26) 


the  last  term  on  the  right-hand  side  of  eq.  (26)  can  be  easily  evaluated  by  forming  a 
pseudo  initial  stress  vector: 


n,k-l 


(27) 


This  pseudo  initial  stress  vector  contributes  to  the  element  load  matrix  in  the  usual 
way.  Because  this  quantity  approaches  zero  as  convergence  occurs,  it  is  not  necessary 
to  evaluate  the  stress-strain  properties  at  the  quadrature  points  as  is  the  case  in 
calculating  the  residual  vector  for  the  tangent  stiffness  method  (see  previous 
discussion).^  The  incluaon  of  a  smilar  term  for  the  tangent  stiffness  method  and  for 
successive  approximations,  evaluated  by  means  of  eq.  (14),  appears  not  to  be  standard 
practice,  and  was  not  considered  in  this  study. 


Althopgh  the  evaluation  of  the  properties  at  the  quadrature  points  is  not  necessary 
to  assure  convergence,  it  might  improve  the  rate  of  convergence;  this  possibility 
was  not  explored  in  the  study. 


Currently,  one  of  the  most  popular  means  for  reducing  the  number  of 
triangularizations  is  to  use  a  quasi>Newton  update  of  the  triangularized  or  inverted 
stiffness  (approximate  3acobian)  matrix.  As  noted  previously,  the  only  quasi^Newton 
methods  considered  in  this  study  are  those  that  directly  update  the  triangularized  form 
of  the  matrix  and  do  not  disturb  its  banded  nature.  In  order  to  satisfy  the  requirement 
of  ease  of  implementatiori,  a  simpler  quau-Newton  update  formula  than  those  available 
in  the  literature  was  sought.  Because  it  is  generally  agreed  that  the  BFGS  update 
formula  is  the  best  available  (see  Fletcher  19S0),  a  formula  of  similar  form  was  desired. 
Denoting  the  triangularized  form  of  the  tangent  stiffness  matrix  as  [Kl*,  the  vqsdate 
for  the  "k"  iteration  expressed  in  index  form  is: 


=  Kf. 
‘Jk-l 


♦  V? 


b, 


I 

£ 

l=l 


«,K*  . 

1  i) 


k-l 


(28) 


where 


a.  s 


b. 


“'i 


1=1 


1=1 


‘k-i*^ 


E  I  K*  6,] 
si^'"l=i  H-l  *’> 


E 

m=i 


(29) 


(30) 


I  =  i  bandwidth  (31) 

The  vector  { A}*  is  the  result  of  a  block  formulation  of  the  residual  vector  {ip}^  ^ 

and  its  reduction  up  to  and  including  row  "i".  The  results  of  the  previous  iteration 
are  denoted  by  {6}  =  j  -  {Au}^j^_2*  '*^hen  eqs.  (29)  and  (30)  are  substituted 

into  eq.  (28)  it  is  easily  seen  that  the  update  formula  has  a  form  similar  to  the  BFGS 
formula  (for  the  Sacobian  not  its  inverse).  It  is  also  easy  to  verify  dwt  the  update 
formula  satisfies  the  "quasi-Newton  condition"  (see  Geradin  et  al  1981).  The 
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implementation  of  eq.  (28)  is  strai^tforward,  with  the  updating  of  K*  proceeding 
simultaneously  with  the  reduction  of  the  current  residual  vector. 


3.5.3  Convergence  Acceleration 

The  acceleration  of  convergence  and  enlargement  of  the  radius  of  convergence 
of  an  iterative  scheme  by  means  of  some  type  of  extrapolation  procedure  can  often 
be  very  cost  effective.  The  most  convenient  way  of  expressing  such  an  acceleration 
is  in  terms  of  a  convergence  (iteration,  relaxation)  factor;  i.e.,  an  improved  estimate 
of  the  solution  vector  { Au}^,  for  iteration  k,  is  expressed  in  terms  of  the  estimate 
obtained  from  the  solution  procedure  { Au}|^,  the  previous  estimate  { Au}j^  j  and  a 
convergence  factor  C.  For  component  "i"  this  extrapolation  is  expressed  in  the  form: 


Aut  =  Au.  *C.  lAu.  -  Au.  1  (32) 

‘k  ‘k-1  *k  *k  ‘k-1 

Written  in  this  form,  a  value  of  Ol  leads  to  "over-relaxation,”  C<1  to  "under- 
relaxation”  and  C=1  to  no  extrapolation.  The  simplest  approach  is  to  use  a  constant 
acceleration  factor  C.  =  C  ,  for  ail  iterations  aiKl  ali  components  of  the  solution 

K 

vector,  and  to  require  the  user  to  sigtply  the  value  based  on  past  experience.  A  more 
successful  approach  is  to  use  some  rational  criterion  to  calculate  near  optimum  values 
of  the  factor  (see  Isaacson  and  Keller,  1966).  For  single  degree  of  freedom  problems 
(unknown  x),  the  most  popular  procedure  is  the  Aitken's  7  method,  which  yields 
extrapolated  estimates  for  odd  iterations  beginning  with  the  3rd,  i.e.. 


K  =  \ 


'‘k  -  Vi*' 


Xk"^ 


Vl  *  \.2 


k=3,5,7... 


(33) 


When  this  equation  is  expressed  in  the  form  of  eq.  (32)  it  gives  a  value  for  the 
convergence  factor  of: 


In  order  to  avoid  possible  numerical  problems  when  one  is  either  far  away  from  or 
very  near  to  the  solution,  it  is  desirable  to  place  limits  on  the  value  calculated  by 
eq.  (34).  For  the  purpose  of  this  study,  such  limits  are  expressed  in  the  form 

I 

where  ^  1  and  is  user  supplied. 

Eq.  (34)  applies  to  a  single  variable;  the  question  is  how  to  extend  it  to  a 

multi-degree  of  freedom  problem.  Two  schemes  are  proposed.  In  the  first,  for  a 

given  iteration,  a  constant  value  of  is  used  for  all  components  of  {Au}|^.  The 

value  of  C,.  is  calculated  by  using  the  norms  (N.  =  £|  Au.  |)  of  the  solution  vector 
K  i  ‘k 

in  eq.  (34).  This  procedure  is  based  on  the  assumption  that  the  convergence 
characteristics  of  all  the  components  of  the  vector  are  similar,  frt  the  second  procedure, 
eq.  (34)  is  applied  to  each  component  of  { Au)  to  give  a  different  convergence  factor 
for  each  component,  for  each  iteration. 

Accounting  For  The  Inconsistency  In  The  Incremental  Stress-Strain  Relation 
As  previously  noted,  until  convergence  occurs,  the  estimates  of 
{AcJ  .  ,  used  in  the  calculation  of  [D1  .  ,  do  not  in  fact  satisfy  eq.  (6).  Three 
methods  for  handling  this  inconsistency  are  explored. 

Because  in  fact  the  inconsistency  disappears  as  global  convergence  occurs  (i.e., 
as  j  *  ^®^nk  2^*  alternative  is  to  do  nothing;  i.e.,  to  rely  completely 

on  global  iteration. 

In  the  second  approach,  local  iteration  is  introduced  in  the  calculation  of  16] 
so  as  to  remove  the  inconsistency  (see  Herrmann  and  Taylor  1974).  Using 


<c<(^ 


and  I  is  calculated.  The  values  of  l  and  l  are  then 

used  in  conjinction  with  eq.  (6)  to  calculate  a  new  estimate  of  stress  {Ac^*.  ,  which 
is  in  turn  used  with  { Ae}^  j  to  calculate  a  new  estimate  for  the  incremental  properties 
[D]*.  This  process  is  continued  trtil  convergence  is  achieved  for  {At}*. 
Because  of  the  global  iteration  which  also  tends  to  remove  the  inconsistency,  the 
convergence  limit  on  the  local  iteration  is  taken  to  be  only  1/10  as  restrictive  as  the 
global  requirement.  The  stress  estimate  is  iteratively  modified  (instead  of  the  strain 
estimate)  in  order  to  maintain  a  compatible  global  displacement  field  as  required  by 
the  admissibility  conditions  of  the  finite  element  procedure.^  The  introduction  of  local 
iteration  (for  all  points  where  the  incremental  stress-strain  properties  are  required)  of 
course  substantially  increases,  in  a  given  iteration,  the  number  of  calls  to  subroutine 
CLAY,  presumably  with  a  corresponding  reduction  in  the  number  of  global  iterations. 

In  the  third  approach,  the  inconsistency  was  expressed  in  terms  of  a  pseudo¬ 
residual  stress  vector,  i.e.  =  ^‘^nk-l  "  ^®^nk-l  ^^®^nk-l’ 

incorporated  into  the  system  residual  stress  vector. 

3.6  Finite  Element  Implementation 

Following  the  instructions  given  in  Herrmann  et  al  (1981),  subroutine  CLAY  was 
installed  in  a  standard  two-dimensional,  four  node,  nonlinear  finite  element  program 
(NTD).  The  program  was  extended  beyond  its  original  successive  approximation 
capability  to  iriclude  the  several  options  outlined  in  Section  3,5.  The  program  was 
used  for  the  numerical  study  described  in  the  following  section.  After  the  completion 
of  the  evaluation,  in  order  to  somewhat  simplify  the  code,  two  of  the  less  robust 
features  were  removed  (see  discussion  in  next  section).  A  listing  and  brief  user's 
manual  for  the  code  are  given  in  Appendix  B. 

^  The  actual  consequences  of  modifying  the  strain  estimate  were  not  studied. 


The  original  NTO  code  was  relatively  efficient,  uncluttered  and  well  documented. 
However,  once  the  several  options  described  above  were  included,  these  characteristics 
were  lost.  The  problem  is  that  the  features  described  in  Section  3.5,  when  taken  in 
all  possible  combinations,  lead  to  a  total  of  72  different  solution  strategies}  henc», 
the  flow  through  the  program  has  become  rather  complex.  Thus,  while  the  code  given 
in  the  Appendix  should  be  of  considerable  value  in  possible  future  extensions  of  the 
numerical  study,  it  is  rwt  recommended  for  extensive  production  applications.  Rather, 
it  is  recommended  that,  in  a  future  project,  a  production  code  with  far  fewer  options 
be  developed. 

«.  NUMERICAL  STUDY 

%.l  Scope  Of  Study 

The  purposes  of  the  project  were  to  compare  the  effectiveness  of  various 
TKimerical  strategies  for  implementing,  in  finite  element  analyses,  the  bounding  surface 
plasticity  nfxxlel  for  cohesive  soils,  and  to  investigate  the  numerical  characteristics  of 
the  model  While  theoretical  considerations  can  lead  to  general  statements  concerning 
convergence,  etc.,  many  subtle  differences  can  only  be  determined  by  numerical 
experimentation.  The  theoretical  foundations  of  the  methods  being  investigated  are 
for  the  most  part  well  established,  and  thus  this  study  concentrated  on  numerical 
experimentation.  Three  representative  problems  were  selected  for  study  arxl  analyzed 
by  the  several  different  solution  strategies  described  in  the  previous  section.  Results 
of  the  several  rmalyses  were  compared  on  the  bases  of  reliability  and  computation 
efficiency. 

4.2  Selection  Of  Representative  Problems 

The  criteria  for  selecting  the  problems  were  that  they  a)  be  related  to  actual 
geotechnical  engineering  problems,  b)  be  numerically  challenging  and  c)  be  simple 
enough  that  the  overaU  cost  of  the  study  would  not  exceed  the  available  resources. 


The  soil  properties  (model  parameters)  used  in  the  study  were  found  by  applying  the 
calibration  procedure  of  DeNatale  (1982)  to  the  test  results  reported  by  3afroudi  (1982) 
(see  also  Herrmann  et  al  1981);  the  properties  are  given  in  Table  1.  The  degree  of 
initial  overconsolidation  (as  a  function  of  depth  in  the  soil  mass)  was  varied  among 
the  several  analyses.  Convergence  of  the  solution  was  determined  by  requiring  that; 

l\tu  -  Au  I 

i  *k  *k-l  <  ERLMT  (36) 

rTASn 

i  k 

A  value  of  ERLMT  =  .001  was  found  to  be  adequate  for  most  cases. 

The  first  problem,  studied  in  Figure  8,  was  a  simple  triaxial  test  requiring  only 
one  element  for  modeling.  Several  hundred  analyses  were  performed.  All  solution 
strategies  were  used,  a  number  of  different  loading  histories  were  evaluated  and  several 
different  initial  overconsolidation  ratios  were  considered.  Because  of  the  simplicity 
of  the  problem,  only  limited  conclusions  could  be  drawn  from  a  comparison  of  the 
results. 

The  second  problem  consisted  of  the  vertical  loading  of  a  strip  footing  resting 
on  a  clay  deposit.  The  grid  used  in  the  analysis  is  shown  in  Figure  9.  The  extent 
of  the  clay  deposit  was  limited  in  size  in  order  to  keep  the  computer  cost  of  an 
individual  analysis  small.^  Analyses  were  performed  for  normally,  moderately  and 
highly  overconsolidated  soils.  For  normally  and  moderately  overconsolidated  soils,  the 
model  is  numerically  so  well  behaved  that  little  distinction  could  be  made  between 
the  several  solution  strategies.  Thus,  the  bulk  of  the  analyses  (18  in  number)  were 
for  a  highly  overconsoli dated  soil.  However,  even  under  these  conditions,  the  problem 


An  analysis  of  a  footing  on  a  highly  overconsolidated  soil  in  which  the  ioad>de formation 
curve  was  taken  beyond  the  maximum  load  required  from  10-15  increments  and  cost 
on  the  order  of  $5-$10. 
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was  not  sufficiently  challenging  to  reveal  major  differences  among  the  several  solution 
schemes  (with  one  exception  which  is  discussed  in  the  following  section). 

The  third  problem  consisted  of  a  highly  overconsolidated  clay  deport  supported 
by  a  rigid  retaining  wall  that  experienced  rotation  about  its  base  (see  Figure  8),  thus 
leading  to  a  agnificant  reduction  in  confining  pressure  adjacent  to  the  wali.  This 
problem  was  found  to  be  more  challenging  and  served  to  distinguish  among  the  several 
solution  strategies.  A  total  of  33  different  analyses  were  performed  on  the  problem.^ 
The  grid  used  for  the  majority  of  these  analyses  is  shown  in  Figure  9.  The  selection 
of  the  grid,  the  convergence  criterion  and  the  solution  increment  size,  required  to 
achieve  accurate  results  involved  performing  analyses  with  both  coarser  and  finer  grids, 
larger  and  smaller  solution  steps  and  different  convergence  criteria. 

The  coiKlusions  drawn  from  a  comparison  of  the  results  of  the  several  analyses 
of  the  three  representative  problems  are  dbcussed  in  the  next  section. 

#.3  Comparison  Of  Results 

The  only  definite  conclusion  that  could  be  drawn  from  problem  #1  was  that  the 
proposal  to  treat  the  inconsistency  in  the  incremental  stress-strain  law  by  calculation 
of  a  pseudo-residual  stress  was  entirely  unsatisfactory  (convergence  could  not  be 
achieved),  and  was  thus  abandoned  for  the  remainder  of  the  study.  A  postmortem 
investigation  of  the  proposed  method  revealed  no  theoretical  justification  for  its  use, 
and  thus  its  failure  is  not  surprising. 

While  the  second  problem  revealed  some  differences  among  the  several  solution 
strategies,  on  the  whole,  convergence  occurred  so  quickly  that,  with  one  exception,  no 
definitive  conclusions  could  be  drawn.  The  one  conclusion  from  the  second  problem 

^  In  addition,  a  significant  number  of  incorrect  runs  were  made  prior  to  the  detection 
of  the  numerical  integration  problem  for  the  residual  vector  described  in  a  previous 
section. 
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was  that  the  quasi-Newton  scheme  is  not  nearly  as  robust  as  the  modified  Newton 
procedure,  and  at  that  point  in  the  study  it  was  dropped  from  further  consideration. 
At  this  stage  in  the  development  of  quasi-Newton  methods,  it  is  concluded  that  the 
several  advantages  claimed  for  optimization  and  nonlinear  elasticity  problems  do  not 
carry  over  to  plasticity  problems,  at  l^st  as  far  as  the  bounding  surface  model  and 
the  scope  of  this  study  are  concerned.  Other  trends  noted  in  the  second  problem  were 
more  clearly  evident  in  the  third  example  and  are  discussed  in  conjurKrtion  with  that 
study. 

Among  the  analyses  performed  on  the  third  problem,  24  used  the  same  grid 
(Figure  8),  solution  increment  size  and  convergence  criterion  and  thus  can  be  directly 
compared.  These  24  solutions  do  not  completely  exhaust  the  32  solution  strategies 
remaining  within  the  scope  of  the  study  after  eliminating  the  two  components  discussed 
in  the  previous  paragraphs  (the  quaa-Newton  scheme  and  the  pseudo-residual  stress 
representation  of  the  inconsistency  in  the  incremental  stress-strain  law).^  It  is  felt, 
however,  that  these  24  cases  sufficiently  span  the  strategy  space  to  be  representative 
and  permit  definitive  conclusions  to  be  drawn.  The  characteristics  of  these  analyses 
are  described  in  Table  2.  In  the  last  column  of  the  table  a  measure  of  the  cost  of 
each  analysis  is  given.  In  the  calculation  of  this  measure,  the  actual  cost  of  the 
analysis  was  slightly  modified  to  account  for  anticipated  future  savings  due  to  improved 
efficiencies  of  the  equation  solver  and  of  subroutine  CLAY,  and,  in  addition,  the 
pre-and  post-processing  costs  are  not  included.  Because  of  differences  in  coding 
practices,  in  reiative  costs  of  computation  and  storage  for  various  computers,  etc., 
these  figures  contain  a  degree  of  uncertainty,  and  thus  differences  of  less  them  25%-50% 


If  consideration  is  given  to  the  considerable  latitude  that  exists  in  assigning  values 
to  (eq.  (35)),  IREPT  (eq.  (25)),  the  step  size  and  the  ratio  of  the  local  to  global 
convergence  criteria  it  is  seen  that  the  actual  number  of  possibilities  is  really  far 
greater  than  32. 
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probably  are  rx>t  significant.  The  lack  of  convergence,  noted  in  the  table,  for  several 
of  the  analyses  was  due  in  some  cases  6hose  with  a  hi^  computational  cost  measure) 
to  the  limit  set  on  the  maximum  number  of  iterations  in  a  given  increment  and  in 
others  (those  with  an  indicate  low  cost  measure)  to  the  satisfaction  of  the  convergence 
criterion,  eq.  (36),  even  though  convergence  had  not  occurred.  In  either  case,  a  lack 
of  robustness  is  indicated.  In  the  following  paragraphs  the  conclusions  drawn  from 
comparing  the  results  of  the  several  analyses  performed  in  this  study  are  discussed. 

One  proposal  often  made  for  nonlinear  inelasticity  problems  is  to  use  small 
increment  sizes  to  avoid  excessive  (all)  iteration.  The  results  of  this  study  do  not 
support  this  suggestion.  In  Figure  10  the  relative  costs  required  to  reach  a  certain 
point  in  the  solution  history  is  plotted  versus  the  ratio  of  the  number  of  steps  used, 
to  the  minimum  number  required  to  accurately  reach  the  point.  This  plot  clearly 
suggests,  from  a  computational  efficiency  standpoint,  that  one  should  use  the  largest 
step  size  that  will  give  acceptable  accuracy.  (If  the  step  size  is  made  too  large,  the 
numerical  integration  error  in  eq.  (6)  becomes  unacceptable.) 

The  methods  of  successive  approximations  and  tangent  stiffness  exhibit  very 
similar  characteristics  with  no  clear  cut  difference  between  them;  this  conclusion  runs 
contrary  to  the  suggestion  of  Owen  and  Hinton  (1980)  of  the  superiority  of  the  tangent 
stiffness  method.  Both  methods,  if  properly  supplemented  with  other  components  of 
the  solution  strategy,  are  quite  robust,  economical  and  are  sufficiently  accurate  for 
solving  bounding  surface  elasto-plasticity  problems  in  geotechnical  engineemg.  The 
method  of  successive  approximations  has  some  slight,  but  not  major  advantages,  in 
ease  of  implementation. 

The  use  of  local  iteration  considerably  improves  the  reliability  of  the  convergence 
criterion  and  in  most  cases  improves  the  efficiency  of  the  analysis. 


The  cxmsequences  of  using  the  modified- Newton  scheme  and/or  an  acceleration 
factor  are  somewhat  interrelated  and  hard  to  separate.^  When  the  modified  Newton 
scheme  is  not  employed,  the  use  of  a  variable  acceleration  factor  is  of  major  benefit 
for  the  method  of  successive  approximations  and  of  lesser  value  for  the  tangent  stiffness 
method.  When  used  in  conjunction  with  the  modified  Newton  scheme,  a  variable 
acceleration  factor  has  some  value  for  the  tangent  stiffness  method,  but  appears  to 
offer  little  advantage  for  successive  approximations.  Of  the  two  schemes  tested  for 
calculating  a  variable  acceleration  factor,  the  one  that  differs  from  component  to 
component  is  preferable.  The  modified  Newton  method  and  the  calculation  of  a  variable 
convergence  factor  are  both  very  simple  operations  to  implement. 

For  the  grids  used  in  this  study,  the  percentage  of  the  total  solution  cost  spent 
in  evaluating  subroutine  CLAY  ranged  from  about  20%  to  60%.  For  the  larger  grids 
needed  for  production  problems,  it  is  anticipated  that  the  costs  would  be  of  the  order 
of  10%-20%. 

5.  CONCLUSIONS 

Based  on  the  comparisons  made  in  this  study  several  conclusions  are  drawn.  To 
what  extent  these  conclusions  are  generally  applicable  to  very  different  geotechnical 
problems  and  to  other  plasticity  models  is  unknown. 

The  bounding  surface  model  for  cohesive  soils  is  simple  to  implement  in  a 
standard,  nonlinear  finite  element  analysis  using  either  successive  approximations  or 

^  This  interrelationship  is  apparently  due  to  the  fact  that  only  occasionally  updating 
the  stiffness  matrix  has  somewhat  of  a  disriptive  effect  on  the  extrapolation  scheme 
used  in  the  acceleration  of  convergence.  This  suggests  that  results  obtained  from 
iterations  involving  updating  and  those  not  involving  updating  should  not  be  used 
together  in  eq.  (34). 


the  tangent  stiffness  method.  It  is  numerically  well  behaved  and  does  not  lead  to 
prohibitive  computational  costs. 

Both  the  methods  of  successive  approximations  and  tangent  stiffness  can  be 
viewed  as  approximations  to  the  Newton-Raphson  method,  and  there  is  little  to  choose 
between  them.  Either  method  would  appear  suitable  for  a  production  finite  element 
program  for  geotechnical  problems.  Sucoessive  approximations  is  somewhat  easier  to 
implement  and  somewhat  more  efficient;  however,  it  is  slightly  less  robust  than  the 
tangent  stiffness  method. 

The  introduction  of  local  iteration  in  the  calculation  of  the  incremental  stress- 
strain  properties  is  desirable,  arxl  wiU  in  the  future  be  incorporated  into  sii>routine 
CLAY  so  as  not  to  clutter  the  logic  of  the  parent  finite  element  program. 

Finally,  it  is  recommended  that  a  production  program  incorporate  provisions  for 
a  ntodified  Newton  analysis  and  for  a  variable  convergence  factor  based  on  the  use 
of  Aitken's  method  applied  to  each  component  of  the  solution  vector.  Both  of 
these  procedures  are  simple  to  implement. 

It  is  concluded  that  the  effective  use  of  a  quasi-Newton  method  for  elasto- 
plasticity  problems  will  need  to  await  the  development  of  simple  but  robust  updating 
formulas  for  the  banded  ipper  triangular ized  form  of  the  stiffness  matrix. 

&  RECOMMENDED  FUTURE  WORK 

Two  major  components  of  the  overall  project  remain;  i.e.,  the  development  of 
one  or  more  production  finite  element  programs  for  geotechnical  engineering  problems, 
and  a  verification  study  for  actual  field  structures  and/or  large  centrifuge  models. 
The  first  of  the  remaining  components  is  further  discussed  in  the  following  paragraphs. 

Prior  to  any  production  program  development,  sii}routine  CLAY  should  be  recoded 
in  order  to  improve  its  efficiency  and  portability.  It  is  then  recommended  that  two 
production  programs  for  geotechnical  engineering  problems  be  developed.  The  first 
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would  be  two-dimensional  and  applicable  to  either  drained  or  undrained  conditions.  It 
would  contain  simple  pre-processing  routines  for  grid  generation  and  the  description 
cf  initial  stress  and  saturation  states.  The  program  would  be  an  outgrowth  of  the 
program  used  in  the  present  study.  The  grid  generation  scheme  and  the  equation  solver 
would  be  replaced  by  recently  developed  optimized  versions.  The  type  of  element 
would  be  carefully  selected,  in  light  of  accuracy  requirements  for  undrained  conditions. 

The  second  program  would  be  three-dimensional  and,  in  addition  to  being 
applicable  to  ideal  drained  and  undrained  conditions,  would  consider  the  effects  of 
partiad  saturation  and  the  movement  of  pore  water  under  saturated  auid  nonhomogeneous 
stress  state  conditions.  Initially,  a  frontal  equation  solver  would  be  used;  time  permitting 
the  new  iteration  scheme  being  developed  by  Professor  Taylor  at  UCB  would  be 
considered  as  an  al terrtative. 

As  time  permits,  special  features  such  as  bending  elements,  incremental 
excavation  and  construction  options,  more  elaborate  pre  and  post-processing  schemes 
and  complete  dyramic  dimensioning  (as  opposed  to  the  partial  dynamic  dimensioning 
used  in  the  current  two-dimensiortal  program)  would  be  included  in  the  programs. 
Carefully  documented  user  manuals  would  be  prepared  for  both  programs. 
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Table  1:  Results  of  the  Manual  and  Computer  Aided 
Model  Calibration  Proceckires  as  Applied  to 
the  Experimental  Data  of  3airoudi  (1982). 


Parameter^ 

Manual 

Automated 

Solution^ 

Solution^ 

« 

X 

0.130 

0.130 

« 

K 

0.018 

0.018 

M  * 
c 

1.18 

1.18 

« 

M 

e 

0.87 

0.87 

v,G 

0.30 

5900 

Pi 

100 

100 

R 

c 

2.40 

2.509 

R 

e 

2.25 

2.246 

0.04 

0.031 

\ 

0.04 

0.034 

T 

0.10 

0.046 

c 

0.71 

0.453 

s 

1.00 

1.000 

h 

c 

2.00 

0.621 

4.00 

0.855 

* 

m 

c 

0.20 

0.200 

« 

m 

e 

0.20 

0.200 

1  a  *  indicates  that  this  parameter  was  assigned  a  fixed 
value  prior  to  using  the  computer  aided  solution 

2  following  discovery  of  an  error  in  the  reduction  of 
the  iaboratory  results,  the  model  was  recalibrated 
manually,  and  reported  earlier  by  Herrmann  et  a! 
1981b. 

3  the  objective  finction  was  created  with  the  options 
set  at  PUM  =  TUM  s  0.20.  All  data  was  weighted 
equally. 
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APPENDIX  A 


NON-HOMOGENEOUS  MODEL  AND  IN  SITU  TEST  RESULTS 

To  date,  the  bounding  surface  model  has  been  used  primarily  to  predict  the 
behavior  of  soils  under  homogeneous  states  of  stress  and  strain.  In  these  applications 
the  model  has  been  shown  to  be  accurate  and  versatile.  However,  in  order  to  fully 
verify  its  predictive  capabilities  and  practical  usefulness,  it  is  necessary  to  incorporate 
the  material  model  into  a  nonlinear  finite  element  code  and  study  its  characteristics 
when  applied  to  nonhomogeneous  laboratory  or  in  situ  boundary- value  problems. 

The  degree  to  which  this  final  step  in  the  model  verification  process  can  be 
completed  is  limited  by  the  rather  sparse  amount  of  published  experimental  data. 
Ideally,  homogeneous  test  results  should  be  provided  to  permit  the  bounding  surface 
model  to  be  fully  calibrated.  Model  test  or  field  measurements  must  also  be  presented 
to  serve  as  a  measure  of  the  formulation's  predictive  capabilities  for  non-homogeneous 
stress-strain  conditions.  Few  published  studies  meet  both  requirements,  although  it  is 
often  possible  to  acquire  all  necessary  data  by  consulting  several  articles  by  the  same 
research  group.  Some  of  the  reported  experimental  data  that  may  be  used  in  this 
final  stage  of  the  verification  process  is  briefly  discussed  in  the  following  paragraphs. 

Ajaz  and  Parry  (1976)  describe  the  response  of  a  compacted  natural  clay  beam 
subjected  to  a  series  of  laboratory  bending  tests.  The  testing  apparatus  was  designed 
to  permit  a  study  to  be  made  of  the  material's  tensile  strength  and  stress-strain 
response  in  uniform  bending.  The  laboratory  results  are  presented  in  the  form  of  a 
moment-deflection  curve  for  the  center  section,  and  by  a  series  of  strain  contours 
corresponding  to  different  applied  bending  moments.  By  assuming  a  particular  set  of 
equilibrium  and  strain-displacement  relations,  the  authors  are  able  to  portray  their 
results  in  terms  of  the  soil's  stress-strain  response.  Plane  strain  conditions  are  assumed. 
Unfortunately,  none  of  the  necessary  soil  properties  are  provided,  and  there  is  not 


sufficient  information  to  calibrate  the  bounding  surface  model.  Additional  aspects  of 
the  authors's  research  are  presented  in  other  published  articles,  and  it  may  be  possible 
to  acquire  the  required  calibration  information  by  consulting  these  sources. 

Baasubramian,  Sivandran  and  Ho  (1979)  describe  the  results  of  three  separate 
full-scale  in  situ  slope  stability  tests  involving  Bangkok  clay.  One  of  these  embankments 
(specifically,  their  "Embankment  I")  could  be  analyzed  with  a  nonlinear,  two-dimensional, 
plane-strain  finite  element  program.  The  experimental  results  are  presented  in  the 
form  of  force-deformation  and  force-pore  water  pressure  histories  for  various  locations 
at  the  test  site.  Although  not  all  of  the  necessary  material  constants  are  provided, 
Bang)<ok  clay  has  been  extensively  tested  in  the  laboratory,  and  therefore  a  full 
complement  of  material  properties  could  be  readily  acquired  by  consulting  additional 
cited  references. 

Hanzawa  (1979)  describes  a  combined  laboratory  and  in  situ  testing  program 
conducted  with  a  natural  clay.  The  study  is  concerned  strictly  with  ultimate  strength, 
and  an  attempt  is  made  to  identify  the  effects  of  such  quantities  as  consolidation 
history,  strain  rate  and  aging  on  shear  strength.  No  complete  stress-strain  response 
is  provided,  and  no  field  or  model  tests  are  reported. 

Desai  et  al  (19S1)  describe  the  results  of  a  series  of  homogeneous  and  model 
tests  on  an  artificial  soil  made  of  oil,  sand  and  clay.  Conventional  triaxial  compression 
and  extension  test  results  are  provided,  together  with  values  of  the  traditional  critical 
state  parameters  X,  k  and  M;  the  bounding  surface  model  could  be  readily  calibrated 
for  normally  consolidated  conditions.  Two  separate  bearing  capacity  tests  were 
conducted  with  scale  models,  and  the  results  are  report  in  the  form  of  force-displacement 
curves.  Both  tests  could  be  analyzed  with  a  nonlinear,  two-dimensional,  plane-strain 
finite  element  program.  Additional  experimental  results  may  also  be  available  in  other 
published  and  referenced  articles. 
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Andersen  and  Stenhamau-  (1982)  reported  the  results  of  three  in  situ  static  plate 
loading  tests  on  heavily  overconsolidated  Haga  clay.  The  OCR  profile  of  the  natural 
deposit  is  given  and  the  experimental  results  are  presented  in  the  form  of  force- 
displacement  and  force  -  pore  water  pressure  histories.  Either  plane  strain  or 
axisymmetric  conditions  are  assumed,  and  the  tests  could  be  analyzed  with  a  nonlinear, 
two-dimensional  finite  element  code.  The  values  of  the  critical  state  parameters  X, 
K  and  M  are  not  provided,  but  could  be  obtained  by  consulting  other  references. 
Otherwise,  there  is  enough  homogeneous  laboratory  data  to  enable  the  Bounding  Surface 
model  to  be  calibrated  for  the  heavily  overconsoli dated  range  and  compressive  stress 
states. 

Radhakrishnan  and  Reese  (1969)  report  the  results  of  a  laboratory  model  study 
in  which  they  studied  the  response  of  homogeneous  and  two-layered  clay  masses  beneath 
a  loaded  strip  footing.  The  experimental  observations  are  presented  in  the  form  of 
force-deformation  histories  at  various  locations  beneath  the  loaded  footing.  No 
consolidation  or  drained  compression  data  is  provided,  and  the  initial  state  of  the 
material  (e^,  ^^)  is  not  specified.  The  index  properties  of  the  two  natural  clays 
are  tabulated,  and  deviator  stress-axial  strain  relations  from  unconsolidated,  undrained, 
(UU)  triaxial  compression  tests  at  three  different  confining  pressures  are  presented 
from  which  some  of  the  necessary  model  parameters  could  perhaps  be  found.  Both 
model  tests  could  be  analyzed  with  a  nonlinear,  two-dimensional,  plane-strain  finite 
element  formulation.  Additional  material  properties  may  possibly  be  found  in  cited 


references. 


APPENDC  B 


USER'S  MANUAL  FOR  NTD 

Nonlinear  2-D  Stress  Analysis  Program  (Pour  Node  Element) 


L.  R.  Herrmann 

Department  of  Civil  Engineering 
University  of  California,  Davis 
August,  19S2 


PART  1:  INPUT! 


The  required  input  data  is  entered  by  means  of  the  following  sequence  of  cards: 


Al.  Title  Card  (18A»)t 

Any  information  that  is  to  be  printed  as  the  title  of  the  problem. 

A2.  Control  Card  (215); 

The  following  card  is  required  to  define  the  desired  analytical  options: 

Columns 

0  -  plane  stress  analysis 
5  MTYPE  =  1  -  plane  strain  analysis 

2  -  axisymmetric  analysis 

10  IHISBF  =  History  function  number  (corresponding  to  the 

history  function  specifications  of  section  Bl) 
for  the  body  for'*':  terms 

A3.  Nonlinear  Aralysis  Control  Card  (515,  2E10.3):  The  following  card  is  required 

to  specify  the  desired  iteration 
options: 


Columns 


1  -  5 

NONLIN 

1  -  successive  approximations 
[2  -  tangent  stiffness  (Newton^  method) 

6  -  10 

ITMAX 

= 

Maximum  number  of  iterations  permitted  in  any 
single  solution  increment 

11  -  15 

NWAY 

=  1 

fo  -  no  local  iteration 
[l  -  with  local  iteration 

16  -  20 

IRPET 

fo  -  reform  every  iteration 

Ik  -  reform  every  K-th  iteration 

21  -  25 

ITFAC 

4 

0  -  no  acceleration 

1  -  constant  acceleration  factor  =  RELAX 

2  -  variable  norm  acceleration 

>3  -  variable  component  acceleration 

26  -  35 

4 

]  for  ITFAC  =  (  0 

RELAX^  J  1 

J  w 

36-45 

ERMAX 

Convergence  criterion  for  the  displacement 
vector  (by  default,  ERMAX  =  0.01) 

Bl.  A  card  with  a  1  punched  in  cohmn  1  followed  by: 

History  Function  Descriptons:  The  following  cards  are  required  for  each  distinct 

function: 

1st  Card  (IX.  1».  15); 

Columns 

2-5  IH  s  Function  number 

6-10  M  s  Number  of  points  needed  to  define 

the  function 

2nd  CardCs)  (8E10.3); 

As  many  cards  as  needed  to  specify  the  M  pairs  of  values  (F  ,  t  ).  The 
initial  card  should  contain  the  values  F,,  t,,  F,,:,,..^  1^,  t..  S:£sequent 
cards,  if  required  (M  >  should  coniaW  thcT  values  Fj,  ty...,  Fj^,  tj^j. 


■jWr, 


B2.  A  card  with  a  2  punched  in  column  1,  followed  by: 


Material 


ies  Array:  The  following  information  must  be  supplied  for  each 
distinct  material: 


1st  Card  (IX.  15.  2E10.3): 

Columns 


2-5 

10 


MNAT 

mrp 


Material  number 


f  1  -  isotropic  linear-elastic 
=  1  2  -  anisotropic  linear-elastic 

^3  -  bounding  surface  plasticity  model  for 
cohesive  soil 


11-20  -S 

21-30  F^,  J 

2nd  Card  (>£10.3): 

magnitudes  of 

the  body  force  components 

Columns  ITYP  =  1 

ITYP  =  2 

ITYP  =  3 

1-10  E 

*^11 

X 

11-20  V 

*>12 

iC 

21  -  30 

*>13 

“c 

31  -  W 

*>14 

*'c 

41-50 

*>22 

51-60 

*>23 

T 

61  -  70 

*>24 

Pf 

71  -  80 

*>33 

Po 

3rd  Card  (8Ei0.3)  —  required  only  if  ITYP  > 

1: 

Columns  ITYP  =  1 

ITYP  =  2 

ITYP  =  3 

1  -  10 

*>34 

"c 

11  -  20 

*>44 

•'c 

21  -  30 

31  -  40 

G  or  V 

41  -  50 

r 

51  -  60 

Pa 

61  -  70 

n  =  M^/M^ 

71  -  80 

58 

P  = 

I  (3E10.3)  —  required  only  if  ITYP  =  3: 


Ol'liili'-’j 


ITYP  =  1 


ITYP  =  2 


3 


B3.  A  card  ^  a  3  punched  in  column  1,  loltowed  by: 


Initial  Stress  Information:  The  following  information  must  be  supplied  for  each 

initial  stress  state: 


1st  Card  (i: 


6E10.3) 


Columns 


2  - 

5 

ISNO  =  Stress  state  number 

6  - 

15 

^  _  vertical  stress  distribution,^  in  the 

16  > 

25 

•ii 

1  =  form  Oy  =  a^  ♦  a^y 

26  - 

35 

"n 

horizontal  stress  distribution,^  in  the 

36  - 

45 

■>21 

form  ^2^ 

46  - 

55 

'<1 

L  _  pore  water  stress  distribution,^  in  the  form 

56  -  65  CjJ 

2nd  Card(s)  (8E10.3) 

[  ^  u  =  Cj  ♦  Cjy 

As  many  cards  as  needed  to  specify  the  M  pairs  of  values  (t  ,  F  ),  The 
initial  card  should  contain  the  values  t.,  F.,  t,*  F,,  .  .  F.. 

Subsequent  cards,  if  required  (M  >  4)  sh6uld  contain  me  values  t.,  F» 
t  F  '  ^ 


Assuming  plane  conditions.  For  axisymmetry,  the  vertical  coordinate  direction  is 
represented  by  z  instead  of  y,  and  the  distributions  are  then  of  the  form  a  =  a,  *  a,z, 
etc.  ^  ^  ^ 


i»  A  card  with  a  4  punched  in  coiumn  1,  followed  by: 


Node 


Point  Array  (IX,  14,  2EI0.3,  15,  3E10.3):  As  many  cards  as  are  necessary 

to  specify  the  locations  of  aU 
nodes  in  the  system: 


Columns 


2 

-  5 

N 

=  Node  point  number 

6 

-  15 

XN 

=  X  -  coordinate  ^ 

16 

-  25 

YN 

=  y  -  coordinate  ^ 

26 

-  30 

INC 

=  Numbering  increment 

i  quantities 

I  associated 

31 

-  40 

D 

=  Spacing  ratio 

S  with  the 

I  straight  and 

41 

-  50 

XC^ 

Coordinates  of  a  i 

=  point  along  the  1 

f  curv^  line 

'  generation 

31 

-  60 

YC  J 

interior  of  the  J 

circular  arc  t  y 

options 

Assuming  plane  conditions.  For  axisymmetry,  x  r  and  y  2. 


4 


B5.  A  card  with  a  5  punched  in  coliann  1,  followed  by: 

Element  Array  (IX.  1»,  M5):  As  many  cards  as  are  necessary  to  specify  all 

elements  in  the  system:t 

Columns 

2  -  3  The  numbers  of  the  four  node 


1 

6 

11 

16 

-  10 

-  15 

-  20 

points  which  describe  the  quaidrilateral 
or  triangulartt  element  (reading  counter¬ 
clockwise  around  the  element) 

21 

-  25 

MN 

= 

Material  number  (corresponding  to  the 
material  descriptions  of  section  B2) 

1 

26 

-  30 

ISNO 

= 

Initial  stress  state  rtumber  (corresponding 
to  the  stress  state  descriptions  of  section 
B3) 

31 

-  35 

NMIS 

= 

Number  of  additional*^ 
elements  in  the  layer 

36 

-  40 

INC 

* 

Numbering  increment 
for  elements  within 
the  layer 

quantities 
associated 
with  the 
element 

i 

41 

-  45 

NMISP 

• 

Number  of  additional 
layers 

generation 

option 

46 

-  50 

INCP 

= 

Numbering  increment  , 
for  layers  ' 

I 


►*  ' 

«►  - 

% 

_ 

^  The  order  of  these  element  cards  need  bear  no  relation  to  the  actual  location  of 
the  elements  within  the  body. 

tt  For  a  triangular  element  the  forth  node  number  is  set  equal  to  the  first. 
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.  B6.  A  card  with  a  6  punched  in  coliann  followed  byt 

Node  Point  Specification  Array  (IX.  I».  2  (13,  12,  Eia3).  E10.3.  215,  2Eia3)t  As 
many  cards  as  are  necessary  to  specify  the  applied  nodal  displacements 
and  loads: 


Columns 
2-5  N 

6-8  IHj 

10  IFj 

11  -  20  Vj 

21  -  23  IHj 

25  IFj 

26  -  35  Vj 

36-45  e 

46-50  N' 

51  -  55  INC 

56  -  65  P|^ 

66  -  75  Pj^y 


Node  Point  number 

History  function  number  (corresponding 
to  the  history  function  specifications  of 
section  Bl)  for  the  1-coordinate  direction. 

“I  indicates  that  an  >PPlied{J“;^„, 
is  specified  in  the  1-coordinate  direction 

Vaiue  ot  the{j"^^„, 
applied  in  the  1-coordinate  direction 


History  function  number  (corresponding 
to  the  history  function  specifications  of 
section  Bl)  for  the  2-coordinate  direction 

indicates  that  an  appUed^Jj^pJ^^^, 
is  specified  in  the  2-coordinate  direction 

vaiue  el  the{j"^^„, 
applied  in  the  2-coordinate  direction 


Angle  On  degrees)  between  the  x.-axis 
and  x(r)-axis 


Final  node  point  in  the\ 
sequence 

Numbering  increment 
for  node  points  within  ^ 
the  sequence 

Values  of  the  pressures 
applied  at  points  N  and  J 
N  "  respectively  ' 


quantities 

associated 

with  the 

boundary 

condition 

generation 

option 
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B7.  A  card  with  a  7  punched  in  column  1,  followed  by: 


One  card  for  each 
history  segment  into 
which  the  incremental 
analysis  is  to  be  sub¬ 
divided:  t 


Columns 

2-5  NMIS  =  Number  of  solution  (or  time)  increments  into 

which  the  history  segment  is  to  be  subdivided 

6-15  TIME  =  Solution  time  at  the  end  of  the  history  segment 

16-25  D  =  Incrementing  ratio  defining  the  time- 

step  lengths  within  the  history  segment  (by 
default,  D  =  1.0) 

B8.  End  Card  (11):  A  card  with  an  8  punched  in  column  1  to  denote  the  end  of 

the  input  data  for  given  problem. 

Cl.  The  above  sequence  of  cards  A1  B8  should  be  repeated  if  additional  analyses 
are  desired. 


PART  n;  EXPLANATORY  NOTES  REGARDING  THE  INPUT 


Nonlinear  Analysis  Strategy  (section  A3): 


The  choice  among  the  available  solution  strategies  is  selected  by  the  specifications 
on  this  card.  The  reader  is  referred  to  Section  3.5  of  the  main  body  of  the  report 
for  a  description  of  the  available  options. 


History  Fmction  Descriptions  (section  Bl); 

The  histories  of  the  body  forces  and  applied  node  point  displacements  and  loads 
are  specified  by  means  of  "history  functions".  These  history  functions  must  belong  to 
one  of  the  following  three  clcisses: 


i)  IH  <  0:  specifies  identical  incremental  values  which  are  equal  to  the 
specified  force/displacement.  The  incremental  Vcdues  are  taken 
to  be  equal  regardless  of  the  relative  lengths  of  the  time  history 
steps  specified  in  section  B7. 


ii)  IH  =  0:  specifies  a  step-function  history  at  time  tj=0;  that  is,  the 
specified  force/displacement  is  applied  entirely  during  the  first 
solution  increment,  and  no  additional  load/displacement  is  applied 
during  the  remaining  solution  increments. 

iii)  IH  >  0:  specifies  the  particular  history  function  IH  defined  in  section 
Bl.  The  form  of  these  functions  is  illustrated  in  Figure  lA  and 
discussed  below. 


Consider  the  case  of  IH  >  0.  At  time  tj=0  the  function  F(tj)  =  need  not 
necessarily  be  zero.  For  a  step  function  load  (at  tj=0  or  at  any  other  time  t^),  the 
history  segment  must  be  described  as  a  very  steep  ramp  (that  is,  =  small  but 

^  0)  in  section  Bl.  The  solution  segment  must  be  similarly  defined  in  section  B7. 
Within  a  particular  history  segment,  linear  interpolation  is  used  to  identify  the  AF 
which  corresponds  to  the  given  time  increment  At.  For  solution  times  beyond  the 
last  specified  point  t^^,  the  final  history  segment  is  extended  indefinitely.  If  a  value 
V  and  a  history  function  number  IH  >  1  is  specified  in  Section  B6  for  some  given 
external  agent,  then  in  the  solution  interval  At  an  incremental  value  of  the  quantity 
equal  to  VAF  is  applied,  where  AF  (Figure  IB)  corresponds  to  history  function  IH. 


Material  Pr( 


ies  Array  (Section 


opert 

For  general  linear,  anisotropic,  elastic  material  behavior  the  stress-strain  relation 
is  of  the  form; 

I  a]  =  [D1  Icl 

For  plane  stress  or  plane  strain  conditions,  the  stress  and  strain  vectors  are  defined 
by; 

t»|T  =  [0^  Oy 

td'^=  t£^ 

while  for  axisymmetry  (with  respect  to  the  z-axis),  the  vectors  become; 
tol^  ,  o, 

ter  =  Ic,  e^  e,  Y„r 

The  21  material  properties  required  to  define  clay  type  materials  (ITYP  s  3)  are  defined 
in  Herrmann  et  al  (19S0)^. 

Initial  Stress  Information  (Section  B3)t 

In  specifying  the  initial  stress  profiles,  it  is  assumed  that  the  coordinate  system 
is  aligned  in  the  following  as  shown  in  Figure  2B. 


Herrmann,  L.  R.,  Dafalias,  Y.F.  and  OeNatale,  3.S.  (1980),  "Bounding  Surface  Plasticity 
for  Soil  Modelling,"  Department  of  Civil  Engineering,  University  of  California,  Davis, 
Final  Report  to  the  Naval  Construction  Battalion  Center,  Port  Huenente,  CA. 


The  program  incorporates  two  data  generation  routines  to  assist  the  user  in 
defining  the  locations  of  the  system's  node  points.  The  use  of  these  options  can,  for 
example,  enable  one  to  describe  the  nodal  geometry  of  an  arbitrarily  large  grid  with 
as  few  as  five  cards.  Note  that  not  all  numbers  between  1  and  the  maximum  node 
number  need  correspond  to  actual  nodes  in  the  body.  For  example,  the  numbering 
scheme  shown  in  Figure  3B  is  permissible,  and  the  coordinates  of  the  non-existent 
nodes  15  and  21  may  or  may  not  be  specified.  This  feature  facilitates  the  use  of  the 
node  point  and  element  generation  options  defined. 

The  straight  line  or  circular  arc  coordinate  generation  option  may  be  used 
whenever  several  sequential  node  points  lie  along  a  single  straight  line  or  circular  arc. 
If  such  a  situation  exists,  it  is  necessary  only  to  enter  the  coordinates  of  the  initial 
and  final  points  of  the  sequence  (denoted  by  N**  and  N,  respectively),  and  the  values 
of  INC  and  D.  The  constant  INC  represents  the  difference  between  any  two  successive 
node  numbers  in  the  sequence,  and  D  defines  the  ratio  of  the  distances  between  any 
two  adjacent  pairs  of  points. 

If,  for  a  node  N,  INC  ^  0,  intermediate  node  points  are  generated  along  a 
straight  line  (XC  =  YC  =  0)  or  a  circular  arc  (XC  <1  0  and/or  YC  /  0)  between  node 
N  and  the  point  described  on  the  preceeding  node  specification  card  N  That  is,  the 
coordinates  of  the  points  N"  +  INC,  N'  ♦  2*INC,  .  .  .,  N  -  INC  are  each  automatically 
found.  A  circular  arc  is  assumed  to  pass  through  the  end  points  of  the  sequence  N'* 
and  N,  and  some  additional  intermediate  point  having  coordinates  (XC,  YC).  This 
intermediate  point  need  not  necessarily  be  a  node. 

The  end  points  of  the  sequence  may  be  entered  in  either  order.  For  example, 
the  segments  illustrated  in  Figure  4B  could  be  defined  by  specifying  the  nodes  in  either 
the  order  7  22  (INC  =  5)  or  the  order  22  7  (INC  =  -5).  The  spacing  of  the 
intermediate  points  (nodes  12  and  17  in  Figure  4B)  is  controlled  by  the  spacing  ratio 
D.  The  segments  shown  in  Figure  4B  could  be  generated  by  specifying  either  the 
order  7  -►  22  and  D  =  2.0  or  the  order  22  7  and  D  =  0.5.  A  value  of  D  =  1.0  would 

result  in  equally  spaced  nodes. 
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The  interior  node  point  generation  option  locates  all  interior  nodes  whose 
coordinates  have  not  been  established  through  the  options  cited  above  (that  is,  all 
points  still  left  undefined  after  the  section  B4  input  has  been  completed).  The  locations 
of  these  undefined  interior  nodes  are  compMJted  by  means  of  the  Taplacian  - 
Isoparametric"  grid  generation  scheme  developed  by  Herrmann  (1976).  In  this  scheme, 
the  coordinates  of  an  interior  node  are  selected  so  as  to  represent  a  weighted  average 
of  the  coordinates  of  the  neighboring  nodes. 

Figure  5B  illustrates  two  grids  that  have  been  prepared  with  the  aid  of  the 
Laplacian  -  Isoparametric  grid  generation  scheme.  Grid  1  was  developed  by  using  the 
straight  line  generation  routine  to  specify  only  the  exterior  (or  boundary)  nodes  (and 
thus,  only  five  cards  were  needed  in  section  B4).  Grid  2  was  developed  in  an  identical 
manner,  except  that  the  straight  line  ger»eration  option  was  also  used  to  define  the 
nodes  lying  along  the  line  3  21.  Note  that  the  exterior  (or  boundary)  nodes  must 

always  be  directly  or  indirectly  specified  in  section  B4. 

Element  Array  (Section  B5); 

If  the  body  can  be  divided  into  layers  of  elements,  and  if  the  material  number 
MN  and  the  initial  stress  state  number  ISNO  is  the  same  for  several  elements  within 
a  iayer  (or,  perhaps,  for  several  layers),  the  node  numbers  of  these  elements  can  be 
established  by  means  of  the  element  data  generation  option.  To  generate  a  sequence 
of  elements  within  a  single  layer,  node  points  are  specified  for  the  first  element  only, 
together  with  appropriate  values  of  NMIS  and  INC  (see  section  B5). 

For  example,  the  bottom  row  of  elements  in  the  grid  of  Figure  6A  could  be 
established  by  entering  either  the  node  numbers  of  element  "A"  and  the  values  NMIS  =  6 
and  INC  =  4  or  the  node  numbers  of  element  "B"  and  the  values  NMIS  =  6  and  INC  =  -4. 
Similarly,  the  left-most  column  of  elements  could  be  established  by  entering  the  node 
numbers  of  element  "A"  and  the  values  NMIS  =  2  and  INC  =  1.  Note,  again,  that  the 
generated  elements  must  be  of  the  same  material  as  the  specified  or»e. 


t  Hermann,  L.R.  (1976),  "Laplacian  -  Isoparametric  Grid  Generation  Scheme,"  Journal 
of  the  Engineering  Mechanics  Division  ASCE,  v.  102,  no.  EM5. 


If  several  layers  of  elements  are  of  the  same  material,  it  becomes  possible  to 
carry  this  option  one  step  further.  For  example,  the  bottom  two  layers  of  elements 
in  the  grid  of  FigMre  6B  could  be  established  by  entering  the  node  numbers  of  element 
"A"  and  the  values: 

NMIS  =  6 

me  =4 

NMISP  =  1 

INCP  =  1 

or,  alternatively, 

NMIS  =  1 

INC  =  1 

NMISP  =  6 

INCP  =  I 

Furthermore,  if  Grid  1  in  Figure  5B  represented  a  homogeneous  body,  the  entire 
element  array  could  be  established  by  means  of  a  single  card  containing,  for  example, 
the  node  numbers  of  element  "A"  and  the  values: 

NMIS  =  3 

INC  =  I 

NMISP  =  4 

INCP  =  5 

Hence,  under  "ideal"  conditions,  the  element  array  for  a  homogeneous  body  could  be 
defined  with  only  a  single  card  in  section  B3. 

Node  Point  Specification  Array  (Section  B6)i 


Boundary  or  interior  node  point  displacement  and  load  specifications  may  be 
given  in  terms  of  either  the  x-y  coordinate  system  (when  0  =  0  in  section  B6)  or  a 
local  Xj-Xj  coordinate  system  (when  6  ^  0),  as  shown  in  Figure  7B.  If  0  =  0,  the 
subscripts  1  and  2  in  section  B6  refer  to  x(r)  and  y(z)  (and  thus  IF^  =  IF^,  etc.)  and 
if  0  ^  0  they  represent  Xj  and  X2  (and  thus  IFj  =  IF^^  ,  etc.). 
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For  each  of  the  two  coordinate  directions,  one  may  specify  either  a  displacement 
(IF  =  1)  or  a  load  (IF  =  0)  by  setting  V  equal  to  the  applied  quantity.  Specified 
displacements  and  loads  are  considered  to  be  positive  when  they  have  the  same  sense 
as  the  positive  coordinate  directions.  If  a  node  is  neither  constrained  nor  subject  to 
an  applied  load  it  need  not  —  and,  for  economy,  should  not  —  be  included  in  the  node 
point  specification  array  (section  B6). 

A  uniform  or  iinearly  varying  pressure  may  be  specified  aiong  a  straight  or 
curved  boundary  (or  an  interior  line)  by  means  of  the  node  point  specification  generation 
option.  To  use  this  option,  the  quantities  IFj,  IF2,  V|,  V2  and  6  in  section  B6  are 
set  to  zero  (or  left  blank),  and  the  appropriate  values  of  N,  N  Pn  and  are 
entered.  For  example,  to  specify  the  boundary  loading  shown  in  Figure  7,  the  user 
would  enter,  on  a  single  input  card,  the  values: 


N" 

11 

N 

2 

INC 

= 

-3 

Pn 

100.0 

Pn- 

s: 

50.0 

Note  that  the  points  N  and  N'*  must  be  specified  in  a  counter-clockwise  order  if  they 
lie  on  an  exterior  boundary  and  in  a  clockwise  order  if  they  lie  along  an  interior 
boundary  (or  "hole").  Note  also  that  the  pressure  specification  cards  must  proceed  all 
other  node  point  specifications  in  section  B6. 

General  Comments: 

1.  It  is  the  responsibility  of  the  user  to  maintain  consistent  units.  The  units  used 
to  describe  the  material  properties  (section  B2)  must  be  consistent  with  those 
used  to  describe  the  initial  states  of  stress  (section  B3),  the  geometry  of  the 
body  (section  B4),  and  the  applied  loading  (section  B6).  The  solution  will  be 
expressed  in  terms  of  units  which  are  consistent  with  those  of  the  input 
specifications. 

2.  Because  the  bandwidth  of  the  simultaneous  equations  is  determined  by  the 
numbering  of  the  nodes,  an  optimal  node  numbering  scheme  is  required  to 
minimize  the  computational  costs  of  a  given  finite  element  analysis.  The 
bandwidth  resulting  from  a  given  numbering  scheme  may  be  computed  in  the 
following  manner: 
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i) 


ii) 

iii) 


denote  the  bandwidth  span  for  any  two  nodes  of  a  given  triangular 
(three  node)  or  quadrilateral  (four  node)  element  as  N.,  where  N- 
is  equal  to  2  plus  twice  the  absolute  difference  in  node  numbers, 
denote  the  maximum  value  of  N|  for  a  given  element  j  as  NE^. 
considering  all  elements  in  the  system,  denote  the  meucimum  value 
NE.  as 


Since  NE  represents  the  bandwidth  of  the  simultaneous  equations,  in 
numbering  the  nodes  it  is  this  quantity  that  should  be  minimized. 


As  the  program  is  now  dimensioned,  the  value  of  must  not  exceed  64, 

the  maximum  node  number  (NPT)  must  not  exceed  900,  the  number  of  elements 
(NELEM)  must  not  exceed  841,  the  number  of  node  point  specifications  (NBPTC) 
must  not  exceed  120,  and  the  number  of  different  materials  (MNAT)  must  not 
exceed  5.  When  changing  the  dimensions  of  the  program  three  separate  areas 
must  be  considered: 

i)  the  COMMON  blocks; 

ii)  the  values  of  KK  =  -f  1  and  LONG  at  the  beginning  of  the 
code;  and, 

iii)  the  dimension  checks  at  the  end  of  subroutine  PREP. 

The  arrays  used  in  the  program  which  must  be  adjusted  to  accommodate  letrger 
analyses  are  related  to  problem  size  in  the  following  manner: 

i)  X(Nj),  Y(N^)  NQ(Nj  +  l),  DlSPLT(2Nj),  SL(2Nj),  SLP(2Nj),  SLPP(2Nj) 

ii)  N0D(N2,4),  MN0(N2) 

iii)  NODBlNj),  BlV(Nj,  3) 

iv)  Q7(LONG) 

v)  PROP(N^,  21),  FXA(N^),  FYA(N^),  ITYPA(N^) 


where: 


N 

N 

N 

N 

N 


1 

2 

3 

4 

5 


maximum  rx>de  number 

maximum  number  of  elements 

maximum  number  of  node  point  specifications 

NEmax 

maximum  number  of  different  materials 


The  size  of  the  constant  LONG  (and  the  dimension  of  the  Q7  array)  must  be 
large  enough  to  satisfy  the  inequality: 


LONG  > 

iTIoX 


whcrre  represents  the  bandwidth  of  the  simultaneous  equations.  Although 

the  value  of  LONG  must  satisfy  the  above  relation  for  f  =  1.0,  it  is  recommended 
to  specify  LONG  on  the  basis  of  f  >  2.0. 
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Axisymmetry  Plane  Strain 


Figure  2B:  Coordinate  Alignments 


LISTING 


fHESET  FREE 

FILE  1 (KIND-DISK,  MAXRECSI7.E=8200,  BLOCESIZE= 1000,  AREAS-20,  AREASIZE=1) 
FILE  2(KIKD-DISK,  rAXRECSIZE^ 150,  BLOCKSIZE- 150,  AREAS=21C,  AREASIZE^I) 
FILE  5{KIND=-READEP,FILETYPE=7) 

FILE  6 (KIND^ PRINTER) 

C 

C  TWO-DIMENSIONAL  NONLINEAR  FINITE  ELEMENT  ANALYSIS 

C 

COMMON  07(8200) 

COIT.'OH  /BLKO/  NPT ,  NELEM ,  NBPTC ,  MTYPE ,  IHISBF ,  ITMAX ,  RELAX ,  ERf ;AX , 

«  VITFAC ,MODNEW, ITFAC , NONLIN , NWAY , IREPET , ITNO , NSTP 

•  ,ITIM,TIMF, TIME, XKP, XXX 
COr3‘ON/BLK2/  X(900), 7(900), NQ(P01 ),DISPLT(l800) 

C0r3 '0N/ELK3/  MNO  (811 1 ) ,  MOD  (8P 1 ,  ) 

CO! 310N /ELKJ< /  NODB (120),BIV(120,3) 

C0n:0N/BI.K5/  SL ( 1800), SLP(  1800 ),SLPP(  1800) 

COMMON /nLE6 /  ROA  ( 4 ) , SCA  (»J ) ,  ETA  ( 4 ) 
con :0N /CNT/  ICNT 1 , ICKT2, ICNT3 

DATA  ROA.SCA.ETA  /I . 0, -1. 0, 1 . 0, - 1 . 0, 1 . 0, -1 . 0, - 1 . 0, 1 . 0, 

»  1.0, 1.0, -1.0, -1.0/ 

c 

C  FORMAT  STATEMENTS 

C 

804  FORMAT (1 1, 14, 2E 1C. 3) 

90C  FORMATdX, 13,  IX, -ITERATIONS  WERE  PERFORMED  AT  TIME',F6.3, 

«  £X,’THE  ERROR  AT  THE  END  OF  THE  PROCESS  WAS',F7.4) 

901  FOPMAT(/, IX, '*»«««  THE  BANDWIDTH  OF’, 17, 

*  2X,  IS  TOO  LARGE  FOR  THE  DIICENSION  OF  --LONG—  ••**•') 

90r  FORMATdX, '«‘^«<^'GONVERGENCE  DID  NOT  OCCUR»*««  ) 

903  FORI  ATdHO) 

C 

C  THE  SIZE  OF  THE  EQUATION  BLOCK  IS  SET 

C 

LONG=8200 

XXX=1.0E+20 

KK-232 

C 

C  INPUT  DATA  IS  READ 
C 

DO  90  I^I.Ki: 

90  ::c(i)^o 

ICNT1-C 
irNT2=0 
:ci.-T3--0 
SLNRt:i.  0.0 
100  IX- 1 

*NOTE***** 

WRITE (6,9000 )ICNT 1 , ICNT2, ICNT5 

9000  FORMAT dX,  NO  TOTAL  SOL=' ,13, IX, 'NO  PARTIAL  SOl= ' ,13, IX, 'NO  CALL 

1  CLAY= '  ,no) 

CALL  PREF(IX,C7(1)) 

C**cf#*n«iJOTE«»*«*«»NOTE««««»«**«»*'‘'« 

ICNTUO 


ooo  of')o  ono  ooo 


ICNT2-0 

1CNT3=0 

IFdX  .EQ.  1)  GO  TO  110 
105  READ(5,604)NSEC 

IF (NSEC  .NE.  0)GO  TO  100 
GO  TO  105 

THE  EQUATION  POSITIONING  MATRIX  FOR  THE  SYSTEM  MATRIX  IS  FORMED 
110  N^2 

DO  160  1=1, NPT 
160  NQ(I-<  1)  =  0 

DO  196  M=1,NELni 
DO  198  J=1,1» 

T=N0D(t',J)n  1 
198  NQ(I)=N 
HQ(1)=1 

DO  200  1=1, NPT 
200  HC(Ih  1)=NQ(I)-:NQ(I-^1) 

THE  BANDWIDTH  -NCOL-  IS  COMPUTED 

NCOL=0 

DO  218  I  =  1,NELEI-1 
DO  202  J=1,A 
JJ=NOD(I,J) 

DO  20?  K=1,4 
Ki:--N0D(I,K)  -  1 
IViNQ(KK)  -  NQ(JJ) 

IF  (NCOL  .LT.  IV)  NCOLsIV 
20?  CO!!TIHUE 
218  CONTINUE 

COMPUTE  MATRIX  SPECIFICATIONS 

HROVJ=NQ(NPT  +  1)  -  1 
L2.-  NCOL- 1 

LUCLONC  -  L2«KC0L)/NC0L 
IF(L1  .GT.  0)  GO  TO  220 
WRITE (6, 901)  NCOL 
GO  TO  105 

2?0  ID7SE= (NPOU-I )/L1+1 

IF(NP0V.'  .CT.  LONC/NCOL)  GO  TO  222 

IDISK=0 

L1:NnOU 

L2=0 

222  LT  LI  L? 

INITIALIZATION 

DO  225  I=1,NR0W 
DISPLT(I)-0.0 
SLP(I)=0.0 
225  SL(I)-0.0 
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ooo  oon  ooo  noon  non 


ITi;i=0 

TIHF^O.O 

MARCHING  TIKE 

2k0  READ(5,8Ci»)  NSEC,NMIS,TIMFS,D 
IF  (NSEC  .!!C.  0)  CO  TO  100 
IF(D  .EQ.  0.0)  D=1.0 

Dl'2-II!:iS 
DUU1.0/DU2 

IF(D  .EQ.  1.0)  CO  TO 

DU1  =  (1.0  -  D)/(1.0  -  D••N^^IS) 

DT=  (TIMES  -  TIHF)»DU1 
DO  760  NSTPS=1,N11IS 
NSTP=f.’STPS 
WPITE(6,903) 

CHAI:GE  DISPLACEMENT  ESTIMATE  AT  BEGINING  OF  NEW  SOLUTION 
SEGirnJT  Ii:  CASE  OF  UNSTABLE  BEHAVIOR 

IF(NSTPS  .GT.  1)G0  TO  255 
DO  250  Ii1,NI?0V.' 

250  SL(I)=-.01»SL(I) 

255  ITir(-ITIM  ^  1 
TIME=TTMF 
TIMF-TI!!F  s  DT 
DT-DT^D 


ITERATION 


ITKO--1 

260  IT!;0--.ITN0  +  1 


DECIDE  IF  SHOULD  UPDATE  STIFFNESS  OR  NOT 


MODr.’EU:- 1 

I- (ITNO-1 )/IREPET 

IFCLSTPS  .EQ.  1  .AND.  ITNO  ,EQ.  0)MODHEW=0 
IFdTNO  .EQ.  I'lREPET  1)i:OBNEl.’  0 

INITIALIZE  SYSTFi:  MATRIX 

DO  270  I=1,nROU 
SLPP(I)-SLP(I) 

SLP(I)=SL(I) 

2',0  SL'I}--0.0 

IFtMODNEV;  .EQ.  1}  GO  TO  300 
LLi:COL«LT 
DO  280  I-1,LL 
ZCO  Q7  :i)--o.c 
300  CONTINUE 
C»"'-f»ctff*t»fjoTE»»««««»«^«'«’NOTE«‘’**««»«» 

IF  :10DNF'..’  .EQ.  OICNTI-ICNTI-.  1 
IFClODllEV;  .EQ.  1  )ICHT2dCHT2-^1 
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non  n  n  o 


CALL  SOLVE  (LI ,L2,LT,HR0W,NC0L, IDISK,Q7,Q7) 

BOUNDARY  POINTS  TRANSFORKiED  TO  X-Y  COORDINATE  SYSTEll 

DO  570  K-I.NBPTC 
KK=IAES(riODE(K); 

N0DB(K)=KK 
AHG--BIV(K,3) 

IF(ANG  .EQ.  0.0)  GO  TO  570 
KUKK/ 1000000 
JJ=NQ(K1) 

IF(JJ  .EQ.  MQ(KUI))  GO  TO  570 
CT--COS(AIIG) 

SA=SIN(ANG) 

D1=£L(JJ) 

D2^SL(JJ- 1) 

SL(JJ)--D1*CT-D2*SA 
SL(JJ-«1)=D1*SA-D2*CT 
570  CONTINUE 

NEED  FOR  ADDITIONAL  ITERATION  IS  CHECKED 

ERNRM. 0.0 
SOLNRMrO.O 
ITSTOP-0 
DO  580  I=1,NR0U 

C«f  f  ♦■*»«§  «note**********note**»****»»*»*»»* 

IF(N0I:LI1?  .EQ.  2)SL(I)--SL(I)-iSLP(I) 

ERNRI5-ABS(SL(I)  -  SLP(l))  +  ERNRF 
SOLNRM=SOLNEt!  -  ABS(SL(I)) 

580  CONTINUE 

C«rr*»«rffrfrijOTE»*»««*»»»»«NOTE«'»»«»« 

IFdTFAC  .NE.  2)G0  TO  58tt 
RELAX:  1.0 

IFCITNO  .EQ.  0  .OR.  (ITNO/2)*2  .HE.  ITNO)GO  TO  582 
CALL  ACCEL (SLNRI'2,SLNRM1  .SOLNRtJ, RELAX, VITFAC ) 

582  SLNRM2:SLNRH1 
SLNRH1=SOLNRM 

584  IF(SOLriRJ1  .LE.  0.0)  SOLMRI1: 1.  OE-20 
ERNRMrERNRIVSOLKRM 
IFdTNO  .LT.  ITI:AX)  GO  TO  590 
WRITE (6, 902) 

GO  TC  600 

590  IF(ERNnil  .GT.  ERt'AX)  GO  TO  650 
600  WRITE(6,900)  ITNO.TIflF, ERNRM 
ITSTOP= 1 
GO  TO  700 

650  DO  670  I:1,NROW 

C*r  *1*1  rrrvf}QYE*cc  If  rt»r 

IFdTFAC  .NE.  3)GO  TO  670 
RELAX  1.0 

IFdTNO  .EQ.  0  .OR.  dTNO/2)«2  .HE.  ITNO)GO  TO  670 
DU-SL(I) 

DUUSLPd) 
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DU2=SLPP(I) 

CALL  ACCEL ( DU  2 , DU  1 , DU , RELAX , VITFAC ) 

670  SL(I)--SLP(I)  4  RELAX* (SL (I)  -  SLP(l)) 

C 

C  STRESS  AND  STRAINS  COMPUTED 

C 

700  CALL  STRESS  (ITSTOP) 

IFdTSTOP  .EQ.  0)  GO  TO  260 
760  CONTINUE 
GO  TO  240 
END 

SURROUTINE  PREP  (IX, SL) 

C 

C  THIS  SUBROUTINE  SETS  OP  THE  DESCRIPTION  OF  THE  PROELEI! 

C 

com  ;0N /BLK  0/  N  PT , NELDl , NBPTC , MTYPE , IHISBF , ITK AX , RELAX , ERMAX , 

•  VITFAC , HODNEW, ITFAC , HONLIN , NWAY , IREPET , ITNO , NSTP 

•  ,  ITTM,  TIflF,  TIME,  XKP,  XXX 

COIfllON /BLK2/  X (900 ) , Y ( 900 ) , HQ ( 90 1 ) , DISPLT (1800) 

C0I2  ;0H /ELF  3/  MHO  ( 84 1) ,  NOD  ( 84 1 , 4  ) 

COmON  /ELK4  /  NODE  ( 1 20 ) ,  BI V  ( 1 20 ,  3 ) 

COM!  :0N  /ELK 7  /  FUN  ( 1 0 , 3  ) ,  FONT  ( 1 0 ,  3 ) ,  N PTS  ( 3 ,  3 ) 

DII^:NSI0N  ANI(12),AN(3),TITLE(20),H0DS(4),IS(8),HN(2), 

•  IIFLG(2),BIVD(2),IH(2),SC0EF(8,6),SL(1),SC0EFS(6) 

C 

C  FORI'JVT  STATEMENTS 

C 

800  F0RMAT(9I5.2E10.2) 

801  F0R.MAT(CE10,3) 

802  FORMAT (I 1,14,915) 

803  FOR!iAT(I1,l4,6E10.3) 

801  FORMAT  (2CA4) 

806  F0RMAT(n,l4,2(I3,I2,E10.3),E10.3,2l5,2E10.3) 

608  F0KMAT(I1,I4,2E10.3,I5,3E10.3) 

900  FORMAT  (1H0  6X  20A4  ///) 

901  FORMAT ( 3 OX, 'PLANE  STRESS  ANALYSIS’,///) 

902  FORMAT  (1H0, //////,  35X, ’•»■•«*  ELEMENT  INFORMATION  •*•••',///, 

1  IX,- ELEMENT ' , 9X , ' ELEMENT ' , 1 5X , ' ELEMENT ' , 

2  12X, 'MATERIAL ',5X, 'INITIAL  STRESS  STATE',/, 

3  IX, 'HUMBER', 10X, 'CENTER', 13X, 'NODE  POINTS', 11X, 

4  'WUnEER',4X,  SIG-V,  6X, 'SIG-H ' ,  8X, 'U '  ) 

903  F0RMAT(30X,  PLANE  STRAIN  ANALYSIS',///) 

904  FORMAT ( 2X, ’SUCESSIVE  SUBSTITUTION  USED  FOP  NONLINEAR  ANALYSIS') 

908  P0P*;AT(2X, 'TANGENT  STIFF  METHOD  USED  FOR  NONLINEAR  ANALYSIS') 

906  FORMAT(I7,6X,2(a6,1PE10.?,a6,I3,8X),3X,A6,0PE10.3,/) 

9C7  FORMAT(30X, ■ AXTSyMf^:TRIC  ANALYSIS',///) 

906  FOnMAT(1HO,9X, ••t**»GEOMETRy»**«*' ,/, 'ONODAL  POINT', 6X, 

•  X-Y  OR  R-r.  COORDINATES',/) 

909  FCRMAT(I11,5X, 1P2E10.2) 

910  F0R»;AT(//,  15X, 'HISTORY  FUNCTION  NO.  M3, /,  17X, 'VALUE' ,  12X, 'TIME' ) 

911  FORMAT  (//, IX, 'ERROR-TOO  MANY  MATERIALS' ) 

912  FORMAT  (//, IX, 'ERROR-TOO  MANY  ELEMENTS') 

915  FORMAT  (//, IX, 'ERROR-TOO  MANY  NODE  POINTS') 
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non 


Ql**  FORMATC//,  IX, 'ERROR-TOO  MANY  NODE  POINT  SPECIFICATIONS ' ) 

915  PORMAT(1KO,//////,10X,*»*«»«NODE  POINT  SPECIFICATIONS*#***' , ///, 

•  UX.'NODE',/) 

916  FORMAT (2X, 'THE  ITERATION  PROCESS  HAS  A  LIMIT  OF',I3.3X, 

•  'ITERATIONS  PER  INCREMENT', 

•  /,2X,'AND  A  CONVERGENCE  REQUIREMENT  OF', F8. A) 

917  P0RMAT(I5,5X,1P2E10.3,AI5,I12,0P3E12.3) 

918  FORMAT (2X, '•••••  DATA  ERROR  —  TOO  MUCH  HISTORY  FUNCTION  DATA', 

•  ’  ••••*',//) 

919  FORMAT{12X, 1PE12.3,5X,E12.3) 

920  F0RMAT(//, '  ERROR-DATA  ERROR  IK  ELEMENT', 15) 

922  FORMAT(2X, 'LOCAL  ITERATION  USED  TO  MODIFYTHE  STRESS  INCREMENT') 
DATA  AHI,AN(3),HN  /AHU-X=,  4HP-X=,  AHU-Y=,  AHP-Y=, 

•  AH0-R=,4HP-R=,AHU-Z=,4HP-Zs, 

•  4H0-1  =  ,  AHP-U,  AHU-2=,  AHP-2=, 

•  AHANC-,  AHIHU,  4HIH2:i  / 

923  FORMATdX,  STIFFNESS  MATRIX  UPDATED  EVERY' ,12,  IX, 'ITERATIONS' , 

I  '  BEGINNING  WITH  SECOND' ) 

924  PORMAT(1X,  A  CONSTANT  ACCELERATION  FACTOR= ' , F5. 1 , IX, 'IS  USED') 

925  FORMATdX, 'A  VARIABLE  ACCELERATION  FACTOR  BASED  OK  THE* 

•  '  SOLUTION  NORM  AND  WITH  LIFtITS  OF’ ,F5.2, IX, 'AND' ,F5.2, 

•  IX, 'IS  USED') 

926  FORMATdX, 'A  VARIABLE  ACCELERATION  FACTOR  BASED  ON  THE' 

•  INDIVIDUAL  SOLUTION  COMPONENTS  AND  WITH  LIMITS  OF', 

•  F5.2,1X,'AKD’,F5.2,  IX, 'IS  USED') 

DATA  IS  /I, 2, 3, 4, 1,2, 3, 4/ 

DX=0.0 

DYiO.O 

PNF=0.0 

NMAT-0 

NPT-0 


INPUT  DATA  IS  READ 

READ(5,805,END=700)  TITLE 
WRITE  (6,900)  TITLE 
READ (5, 600)  MTyPE,IHISBF 

READ  (5,800  )r;OHLIN ,  ITMAX ,  NWAY ,  IREPET ,  ITFAC ,  VITFAC ,  ERMAX 
IF(MTYPE  .EQ.  0)  WRITE (6, 901) 

IF(MTYPF  .EQ.  1)  WRITE(6,903) 

IFdtTYPE  .EQ.  2)  WRITE(6,907) 

XKPrO.O 

IF(MTYPF.  .EQ.  0)  XKP=1.0 

IFdTMAX  .LE.  0)  ITMAXil 
TF(ER!:AX  .LE.  0.0)  ERMAX-0.01 
IF(RELAX  .LE.  0.0)  RELAX=1.0 
IFdlONLIN  .EQ.  1  )WRITE(6,904) 

IF(HONLIN  .EQ.  2)WRITE(6,905) 

IF(NUAY  .EQ.  1)WRITE(6,922) 

IFdREPET  .EC.  0)TREPET-1 

HRITE(6,923)IREPET 

WRITE (6, 9 16)  ITMAX, ERMAX 

IF(VITFAC  .EQ.  0. 0)VITFAC- 1. 0 

IFdTFAC  .EQ.O  .OR.  VITFAC  .EQ.  1.0)G0  TO  8 


IFCITFAC  .GT.  1 )G0  TO  5 
RELAX=VITFAC 
WRITE{6,921»)RELAX 
GO  TO  8 

5  DU=1.0/VITFAC 

IFCITFAC  .EQ.  2)WRITE(6,925)VITFAC,DU 
IFCITFAC  .EQ.  3)WRITEC6,926)VITFAC,DU 
8  READC5,80r)  NSEC 

INPUT  CONTROL  UNIT 

10  GO  TO  C20,60,75,85, 112, 175, 195), NSEC 

TIliE  FUNCTIONS 

20  READC5,80?)  nsec, I, K 

IFCNSEC  .HE.  O)  GO  TO  10 

IFCI  .LT.  .AND.  N  .LT.  11)  GO  TO  22 

V?RITEC6,918) 

GO  TO  TOO 
22  NPTSCI, 1)=C 
NPTSCI,2)=1 
KPTSCI,3)^N-1 

READC5,8C1)  CCFUHCJ,I),FUHTCJ,I)),J=1,N) 
WRITEC6,910)  I 

WRITE C6, 9 19)  CCFUNCJ,I),FUNT{J,I)),J=1,N) 
GO  TO  20 

MATERIAL  PROPERTIES 

60  ItUI 

CALL  PROPTyClN,NZ,NMAT,HSEC) 

GO  TO  10 

INITIAL  STRESS  SPECIFICATIONS 

75  READC5,803)  NSEC,I, CSC0EFSCJ),J=1,6) 
IFCNSEC  .NE.  0)  GO  TO  10 
DO  77  J=1,6 

77  SCOEFCI, J)-SCOEFSCJ) 

GO  TO  75 

NODE  POINT  COORDINATES 

85  READC5,806)  NSEC,N,XP, YP,INCR,D,XC,YC 
IFCNSEC  .NE.  0)  GO  TO  10 
XCK)=XP 
YCN)=YP 
KQCN)=-2 

IFCINCR.EQ.C)  GO  TO  111 
IFCD.EQ.0.0)  D=1.0 
NM-CN  -NS)/INCn 
NMIS=IABSCKi:) 

IHCR=INCR«NH/Nt;iS 


D02=W^IS 

DUU1.0/DU2 

IF(D.EQ. 1.0)  GO  TO  87 

DU1=(1.0-D)/(1.0-D»»NHIS) 

87  IF(XC  .EQ.  0.0  .AND.  IC  .BQ.  0.0)  GO  TO  95 

GENERATE  POINTS  ON  ARC 

C1U2.0»(XC-XS) 

C12=2.0»(YC-YS) 

C21=2.0«(XC-XP) 

C22=2.0«(YC-YP) 

DUa1.0/(C11»C22-C12«C21) 

B 1 =XC«XC-XS»XS+YC»YC-YS«YS 
B2=XC«XC-XP»XP+YC»YC-YP»YP 
X0= (C22«B 1 -C 1 2»B2 )«DU 
Y0= (C 11 •B2-C2 1 "B 1) *00 
RC=SQRT( (XC-X0)«»2+(TC-Y0)»«2) 

TH0=ATAN2( (YS-YO), (XS-XO) ) 

IF(THO  .LT.  0.0)TH0=6.2831853+TH0 
DY= - (XP-XO )»SIN (THO )+(YP-YO ) *COS (THO ) 

DX=  (XP-XO )»C0S(TH0)4 (YP-Y0)«SIN (THO) 

DTH=  ATAM2(DY,DX) 

ZCs (XC-XS  )• (YP-YS )- (XP-XS )• (YC-YS ) 

IF(ZC  .GT.  0.0  .AND.  DTH  .LT.  0.0)DTHs6.283l853+DTH 

IF(ZC  .LT.  0.0  .AND.  DTH  .CT.  0. 0)DTH=-6. 2831653+DTH 

DTHiDTH»DU1 

DO  90  I=2,NHIS 

THO=THO+DTH 

NS=NS+INCR 

NQ(NS)=-2 

X(NS)=XO+RC»COS(THO) 

Y(NS)=YO-^RC»SIN(THO) 

90  DTH=DTH»D 
GO  TO  111 

GENERATE  POINTS  ON  STRAIGHT  LINE 

95  DX=(XP  -XS)»DU1 
DY=(YP  -YS)»DU1 
DO  110  I=2,NMIS 
NS=NS+INCR 
XS=XS+DX 
YS=YS-tDY 
NQ(NS)=-2 
X(NS)--XS 
Y(NS)=YS 
DX=DX»D 

110  DY=DY*D 

111  XS-XP 
YS=YP 
NSsN 

IF(NPT  .LT.  NS)  NPT=NS 
GO  TO  85 


ELEMENT  INFORMATION 


112  N=1 

1 1 3  READ (5,802)  NSEC , ( NODS ( I ) , I = 1 , 4 ) , MWOR , ISIGN , NMISP, INCRP , NMIS , INCR 
IF(NSEC  .NE.  0)  GO  TO  132 

DO  116  1=1, 4 
116  N0D(N,I)=?10DS(I) 

NS=N 

INCRS=0 

INCRZ=INCR 

NMISZ=NMISP 

MN0(N)=MII0R«100  ISIGH 
120  DO  125  K=1,A 
125  N0D(K,?1)=N0D(NS  ,M)+INCRS 
Mf!0(N)=MII0(NS) 

N=N+1 

INCRS=INCRS+INCRP 

HMISP=fJMISP-1 

IF(NMISP  .GE.  0)G0  TO  120 

HH1SP=KMISZ 

INCRS--INCRZ 

INCRZzINCRZflHCR 

NMIS=NMIS-1 

IF (NMIS  .GE.  0)G0  TO  120 
GO  TO  113 

GENERATE  COORDINATES  FOR  UNSPECIFIED  INTERIOR  NODES 

132  NELEM  =N-1 
FAC1T=1.3 

DETERMINE  WHICH  ELEMENT  SURROUND  EACH  NODE  AMD  HOW  MANY  NODES 
NEED  TO  BE  GENERATED  AND  MAKE  STARTING  ESTIMATE  FOR  THEIR  COORDINATES 


KMIS=0 

L0C=1 

DO  139  K=1,NPT 

IF(NQ(K)  .LT.  0)  GO  TO  139 

N0C=0 

DO  1-7  1=1, NELEM 
DO  135  J=1,4 
JJ=J 

135  IF(K  .EQ.  NOD(I,J))GO  TO  136 
GO  TO  137 

136  K=10»I+JJ 
SL(LOC)=FLOAT(N) 

LOC:LOCt 1 
NOC=NOCf 1 

137  CONTINUE 

IFdJOC  .GT.  0)G0  TO  138 

HC(K)=-1 

GO  TO  139 

136  N=(LOC-NOC)«100-NOC 


NMIS=NMIS+6 

MQ(K)=N 

IF(K.LT.3)  GO  TO  139 

X(K)-0.5«(X(K-2)+X(K-1)) 

Y(K)=0.5«(Y(K-2)+T(X-1)) 

139  CONTIMUE 

IF(NMIS.EQ.O)  GO  TO  157 
WTL=0.0 

ITERATE  TO  LOCATE  OWSPECIFIED  NODES 
AS  A  WEIGHTED  AVERAGE  OF  NEIGHBORS 

DO  155  NN=1,NMIS 
IOT=0 

DO  150  J=1,HPT 
N=NQ(J) 

IF(N  .LT.  0)  GO  TO  150 
LOC=N/100 
NOC=MOD(N, 100) 

WT=0.0 

XS=0.0 

YS-0.0 

DO  1A0  JJJ=1,NOC 
N=IFIX(SL(LOC)iO.  1) 

LOC=LOC+1 

IE=N/10 

JJ=K0D(N,10) 

I=IS(JJ43) 

I=K0D(1E,I) 

K-IS(JJ->-1) 

K=NOD{IE,K) 

L=IS(JJ-*-2) 

L=NOD(IE,L) 

XS=XS  ♦X(I)  ♦X(K)  4.X(L)»WTL 

r3=YS  4Y(I)  +Y(K)  4Y(L)»WTL 

1A0  WT=WT4-2.04.WTL 
D1=X(J) 

D2=Y(J) 

X(J )=( 1. 0-FACIT)»D1+FACIT«XS/VT 
Y  ( J )  =  ( 1 . 0-FACIT )  •D24-FACIT*  YS/WT 
D1=ABS((X(J)-D1)/(ABS(D1)4-1.0E-20)) 
D2=ABS((Y(J)-D2)/(ABS(D2)+1.0E-20)) 
IF(Dl4D2  .GT.  .0001)IOT=1 
150  COKTIIiUE 

IFdOT  .EQ.  0)  GO  TO  157 
WTL--1. 

155  CONTINUE 

PRINT  NOD’!  AND  ELEMENT  DATA 
157  WRITE(6,906) 

WRITE(6,909)  (N,X(N),Y(N),N=1,NPT) 

WRITE(6,902) 

DO  168  N=1,NELEM 


ooo  non  oooo 


MIIORS=MH0(N)/100 
ISIGNrMOD (MNO (N ) , 1 00 ) 

MTJ0(N)=KN0HS 

CHECK  FOR  NEGATIVE  ELEMENT  AREA,  INITIALIZE  STRESSES  AND  STRAINS, 
PRINT  INFORMATION,  FORM  ISOPARAMETRIC  TRANSFORMATION  AND  STORE 

N1=N0D(N, 1) 

N2=H0D(N,2) 

N3=NOD(N,3) 

Nl»=N0D(N,'*) 

A1=X(Nl)*(Y(N2)-Y(N»»))-*-X(N2)»{Y(Ni4)-y(Nl))+X(NJ<)*(y(N1)-y(N2)) 
A2=X(H2)«(Y(N3)-Y(Nl«))+X(H3)*(Y(K4)-Y(N2))+X(NH)»(Y(N2)-y(N3)) 
IF(A1+A2  .GT.  0.0)  GO  TO  166 
IX=0 

WRITE(6,920)  N 
GO  TO  167 

166  CALL  CEOt;(N,MNORS,ISIGN,SCOEF,SIGH,SIGV,U,XC,YC) 

167  WRITE(6,917)  N,XC,YC,N1,N2,K3,N4,MN0RS,SIGV,SIGH,U 

168  CONTINUE 
GO  TO  10 

NODE  POINT  SPECIFICATIONS 

175  WRITE  (6,915) 

1=1 

160  READ(5,806)  NSEC,KK,  (IH(N),IIFLG(N),BIVD(N),N=1,2),TH,KKP,INCF., 
•  PJ,PK 

IFCNSEC  .NE.  0)  GO  TO  10 
IADD=0 

if(i:type  .EQ.  2)  IADD=4 
IF(TH  .HE.  O.O)  lADDsB 
NMIS=0 
Nr=1 

IF(KKP  .EQ.  0)  GO  TO  l85 
IFdNCR  .EQ.  0)  INCR=1 
NH=(KKP-KK)/INCn 
m'is=iABs(mi) 

INCR=INCR»NM/NMIS 

N1-U  1  -*  Nf'-IS 

DP=PK-PJ 

M=KK 

XL-- 0.0 


GENERATE  SPECIFICATIONS  FOR  INTERMEDIATE  NODES 


DO  183  L=1,NriIS 
MP=M 


M=MdNCP. 

183  XL=SORT((X(!;)-X(MP)  )«*2*(Y(K)-y(MP)  )»*2)+XL 
DP=DP/XL 


RJ=1.0 

IFCMTYPE  .EQ.  2)  RJ=X(KK) 
185  PXB:.0.0 


89 


PYB=0.0 

DO  190  L=1,NH 

PNBrO.O 

IF(L  .EQ.  MM)  GO  TO  186 
DXa(KK  INCR)  •  X(KK) 

DYsKKK  +  INCR)  -  Y(KK) 

PKsPJ  ♦  DP»(SQRT(DX»DX  +  DY*DY)) 

RK=1.0 

IFdITYPE  .EQ.  2)  RKrX(KK  ♦  INCR) 

PHB=(3.0»PJ»RJ  ♦  PJ»RK  PK»RJ  +  PK»RK)/12.0 
PNF=(PJ«RJ  *  PJ«RK  ♦  PK*RJ  ♦  3.0«p|C«RK)/12.0 

CALCULATE  EQUIVALENT  NODAL  FORCES  DUE  TO  SPECIFIED  PRESSURE 


186  BIV(I,1)=EIVD(1)  +  PXB  -  DY»PNB 
EIV(I,2)=BIVD(2)  +  PYB  +  DX«PNB 
PXB=-DY«PKF 

PYB=  DX»PNF 

PJ=PK 

RJ:;RK 

BIV(I,3)*TH 

NODE(I)=KK» 1000000  ♦  IK(1)»10000  ♦  IH(2)»100 
•  ♦  IIFLG(1)«10  ♦  IIFLG(2) 

DO  187  J=1,2 
K=2*J-UIADD 
AN(J)=ANI(K+1) 

IF(IIFLG(J)  .EQ.  0)  GO  TO  187 
AU(J)=ANI(K) 

187  COllTINUE 

WRITE(6,906)  KK, (AN(J),BIV(I, J),HK(J),IH(J), J=1,2),AN(3),BIV(I, 3) 
KK=KK+IHCR 

BIV(I,3)=TH»0.017‘*533 
NBPTCrl 
190  1=1+1 

GO  TO  180 


THE  SIZE  OF  THE  PROBLEM  IS  CHECKED  TO  SEE  IF  IT  IS  TOO  LARGE 
AND  DATA  ERRORS  ARE  SOUGHT 


195  IF(NELEM  .LT.  842)  GO  TO  206 
WRITE  (6,912) 

IX=0 

206  IFdlPT  .LT.  901)  GO  TO  207 
WRITE  (6,913) 

IX=0 

207  IF(NBPTC  .LT.  121)  GO  TO  206 
WRITE  (6,914) 

IX=0 

206  IF(NHAT  .LT.  4)  GO  TO  210 
WRITE  (6,911) 

IX=0 

210  RETURN 
700  STOP 


ooo  nor>o  ooo 


SUBROUTINE  PROPTY (IN ,MH , NMAT, NSEC ) 

THIS  IS  THE  MASTER  SUBROUTINE  FOR  SUPPLYING  MATERIAL  PROPERTIES 

COWION /ELKO/  NPT , NELEM.NBPTC .MTYPE , IHISBF, ITMAX, RELAX , ERMAX , 

•  VITFAC .MODNEW, ITFAC , NONLIN , NWAY , IREPET, ITNO , NSTP 

•  ,ITIM,TIMF, TIME, XKP, XXX 
COrUlON/BLKI/  PROP(3,21),FXA(3),FYA(3),ITYPA(3) 

C0W:0I1/BLK8/  C(4,4),S0(M),XPV(5),YPV(5),XJC0B(5),C1(8,8),ZY(8), 

•  FV(5,4),GV(5,4),XKV(5,4),SIGT(4),DSIG(4),EPT(4),DEP(4), 

•  STOR(6),PWPT,DPWP,CSC(4,4) 

COMMON/CNT/  ICNT1,ICNT2,ICNT3 

DIMENSION  SIG3D(6),DSIG3D(6),EP5D(6),DEP3D(6),C3D(6,6),CB3D(6,6) 
801  F0RMAT(8E10.3) 

807  F0RMAT(I1,I4,I5,4E10.3) 

904  F0RMAT(//,  IX, 'MATERIALM3.2X, ’IS  ISOTROPIC  WITH  FX  (R  )=  ’  ,E10.  3,4X, 

•  'FY(Z)=',E10.3,4X,*E  =’,E10.3,4X, 'AND  POISSONS  RATIO  =',F5.2,/) 

905  F0RMAT(//, IX, 'MATERIAL', 13, 2X, 'IS  ANISOTROPIC  WITH  FX(R)=' , 1PE10. 3 

•  ,  4X,  'FY(Z)  =  ',1PE10.3,4X,'AND',/,5X,'Cn=',iPE10.3, 

•  5X, 'C12=',0PE10.3,5X, 'C13=' ,0PE1C.3,5X, 'C14= ' , OPEIO. 3, /, 

•  24X, 'C22=',1PE10.3,5X, 'C23=’ ,0PE10.3,5X, 'C24= • , OPEIO. 3, /, 

•  43X,’C33=',1PE10.3,5X,’C34=’,OPE10.3,/, 

•  62X, 'C44=', 1PE10.3) 

GO  TO  (100, 300), IN 

ffffttfoffvtc  read  materials  PR0PERTIES»*»*»»*»*»*« 

100  READ(5,807)  NSEC,I,ITYP,FX,Fy 
IF (NSEC  .NE.  0)  RETURN 
IF (NMAT  .LT.  I)  NMAT=I 
FXA(I)=FX 
FYA(I)=Fy 
ITYPA(I)=ITYP 
GO  TO  (150, l60,200),ITyP 

ISOTROPIC  ELASTIC 

150  READ(5,801)E,GNU 

V’RITE(6,904)  I,FX,FY,E,GHU 
DU=E/((1.0+GNU)  •  (1.0-2.0»GNU)) 

Cl 1:DU«( 1.0-GNU) 

C12=DU*GNU 

C13=C12 

C22=C11 

C23=C12 

C33=C11 

C44=E«0.5/(1.0+GHU) 

C14iC.0 
C24=0.0 
C34-0.0 
GO  TO  165 
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onooo  ooooooooo  oo 


ANISTROPIC  ELASTIC 

160  HEAD(5,801)C11,C12,C13,C1A,C22,C23,C2»»,C33iC3'I,C4>I 

WRITE(6,905)  I,FX,FY,C11,C12,C13,C1A,C22,C23,C2A,C33,C34,CA4 

165  PR0P(I,1)=C11 
PR0P(I,2)=C12 
PR0P(I,3)=C13 
PR0P(I,^)=C1A 
PR0P(I,5)=C22 
PROP(I,6)=C23 
PROP(I,7)=C2»4 
PROP(l,8)=C33 
PROP(I,9)=C3‘» 

PROP  (I,  10)=Cl»l» 

GO  TO  100 

BOUKDING  SURFACE  PLASTICITY  FOR  COHESIVE  SOIL 

200  CALL  RPR0P(PR0P(I,1)) 

00  TO  100 


••#»«««»#»«CALCULATE  INCREr^NTAL  PROPERTIES 

300  ITYP=ITYPA(MH) 

GO  TO  (250,250,A00),ITYP 

LINEAR  ELASTICITY 

250  IF(ITIM+ITNO  .GT.  2)RETORN 
C(1, 1)=PR0P(HN, 1) 

C(l,2)=PROP(MN,2) 

C(1,3)=PROPOnJ,3) 

C(1,4)=PR0P(15N,4) 

C(2,2)=PR0P(I1N,5) 

C(2,3)=PPOP(KN,6) 

C(2,4)=PR0P(MN,7) 

C(3,3)=PROP(MN,8) 

C(3r4)=PROP(MN,9) 

C(4,4)=PR0P(MN, 10) 

S0(1)=0.0 
S0(2)=0.0 
S0(3)=0.0 
S0(4)r0.0 
DO  280  J=1,4 
DO  280  K=1,J 
280  C(J,K)=C(F,J) 

RETURN 

BOUNDING  SURFACE  MODEL  FOR  COHESIVE  SOIL 

400  RTlsl.O 
RT2s1.0 


CHANGE  SIGN  OF  STRAIN  ESTIMATE  AT  BEGINING  OF  NEW  SOLUTION 
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D-ni24  see  nuhericrl  inplehentotion  of  the  cohesive  soil  bounding 

SURFACE  PLflSTICITV.  .  (U)  CflLIFORNIS  UNIV  DflVIS  DEPT  OF 
CIVIL  ENGINEERING  L  R  HERRHflNN  ET  RL.  FEB  83 
UNCLASSIFIED  NCEL-CR-83.  818  Ne2474-82-C-827e  .  F/D  8/13 


ooooo  oonoooooo  oo 


ANISTROPIC  ELASTIC 

160  READ(5,801)C11,C12,C13tClA,C22,C23.C24,C33tC34«C41l 

WR1TE(6,905)  I,FX,FT,C11,C12,C13,C14,C22,C23,C24,C33,C34,C44 
165  PR0P(I,1)=C11 
PR0P(I,2)=C12 
PR0P(I,3)=C13 
PROP(I,4)sCl4 
PR0P(I,5)sC22 
PROP<l,6)=C23 
PROP(I,7)=C24 
PROP(I,8)=C33 
PR0P{I,9)=C34 
PROPd,  10)=C4l| 

GO  TO  100 

BOOMDING  SURFACE  PLASTICITY  FOR  COHESIVE  SOIL 

200  CALL  RPR0P( PROP (1,1)) 

GO  TO  100 

•••••••••••CALCULATE  IHCREUBNTAL  PROPERTIES^^^^^»^^»^*» 

300  ITYP=ITYPA(MM) 

GO  TO  (250,250,400),ITYP 

LINEAR  ELASTICITY 

250  IFdTIM+ITNO  .GT.  2)RETURN 
C(1,1)sPROP(MN,1) 

C(1,2)=PROP(MK,2) 

C{1,3)=PROP(»aJ,3) 

C(1,4)=PROP(WJ,4) 

C(2,2)=PR0P(HN,5) 

C(2,3)=PR0P(r^N,6) 

C(2,4)=PROP(MW,7) 

C(3,3)=PR0P{I3l,8) 

C(3,4)=PR0P(KN,9) 

C(4,4)=PROP(MN,10) 

S0(1)s0.0 
SO(2)sO.O 
S0(3)=0.0 
S0(4)=0.0 
DO  280  Js1,4 
DO  280  Kr1,J 
280  C(J,K)=C(li:,J) 

RETURN 

BOUNDING  SURFACE  MODEL  FOR  COHESIVE  SOIL 

400  RTUI.O 
RT2=1.0 

CHANGE  SIGN  OF  STRAIN  ESTIMATE  AT  BECINING  OF  NEW  SOLUTION 
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C  SEGMENT  IN  CASE  OF  UNSTABLE  BEHAVIOR  AT  END  OF  PREVIOUS  ONE 
C 

IF(NSTP+ITNO  .GT.  t)GO  TO  H05 

RT1=.01 

RT2=-.01 

405  DO  410  1=1,4 

SIG3D(I)=-SIGT(I) 

DSIG3D(I )=-DSIG(I )»RTl 

DSIG(I)=DSIG(I)»RT1 

EP3D(I)=-EPT(I) 

DEP3D(I)=-DEP(I)»RT2 
410  DEP(I)=DEP(I)»RT2 
DPWPT=DPWPT»RT1 
DO  415  1=5,6 
DSIG3D(I)=0.0 
SIG3D(I)=0.0 
EP3D(I)=0.0 
415  DEP3D(I)=0.0 
ITNOP=ITNO 
LITNO=0 

420  LITN0=LITN0+1 
ITN0P=ITN0P+1 

Ct«»ii«*«*»«NOTE«*«*««««***NOTE*»«« 

ICNT3=ICNT3+1 

CALL  CL AT ( 3 , ITIM , ITNOP , PROP (MN , 1 ) , STOR , SIG3D , EP  3D , DSIG  3D , DEP3D 
1  C3D,CB3D,/'i/PT,DPWPT,CAM,1,0) 

R1  =  .5 

IF(NONLIN  .EQ.2)R1=1.0 
R2=1,0-R1 

IFCNOHLIN  .EQ.1  .AND.  MODNEW  .EQ.  1 )G0  TO  460 
DO  450  1=1,4 
S0(l)=0.0 
DO  450  J=1,4 

CSC(I,J)=0,5«(C3D(I,J)+CB3D(I,J)) 

450  C(I,J)=C3D(I,J)»R1+CB3D(I,J)»R2 
460  IF(NWAY  .LT.  1 )G0  TO  625 
DO  560  1=1,4 
DU=DSIG(I) 

DO  550  J=1,4 

550  DU=DU-0.5*(C3D(I,J)+CB3D(I,J))»DEP(J) 

560  SO(I)=DU 
ER=0.0 
XNRM=0.0 
DO  600  1=1,4 
ER=ER+ABS(SO(I)) 

XNRM=XI.’RM+ABS  (DSIG  (I ) ) 

DSIG(I)=DSIG(I)-SO(I) 

600  DS1G3D(I)=-DSIG(I) 

IF(XKRM  .EC.  0.0  )G0  TO  625 
IF(LITNO  .EQ.  ITMAX)GO  TO  620 
IF(ER/XNRM  .GT.  10.»ERMAX)GO  TO  420 
620  DU=ER/XHRH 

625  IF(MODNEW  .EQ.  0  .OR.  NONLIN  .EQ.  2)G0  TO  629 
DO  628  1=1,4 


S0(l)=0.0 
DO  628 

628  SO(I)=SO(IK(0.5»C3D{I,J)+0.5«CB3D(I,J)-C(I,J))»DEP(J) 

RETURN 

C*«t»*«i*«*NOTE**«*«**«**NOTE 

629  R 180.5 

DO  630  l8l,H 
S0(I)=0.0 
DO  630  J8l,4 

630  SO(I)=SO(I)+(R1-0.5)»{C3D(I,J)+CB3D(I,J))»DEP(J) 

RETURN 

END 

SUBROUTINE  STRESS (ITSTOP) 

C 

C  THIS  SUBROUTINE  CALCULATES  AND  PRINTS  ELEMENT  STRESSES  AND  STRAINS 
C 

COI^ON /BLKO/  NPT , NELEM. NBPTC , MTYPE , IHISBP, ITNAX , RELAX , ERMAX , 

•  VITFAC , MODNEW , ITPAC , NONLIN , NVA Y , IREPET , ITNO . NSTP 

•  .ITIM.TIMF, TIME, XKP, XXX 

COMMON /ELK2/  X (900 ) , Y (900 ) , NO (90 1 ) , DISPLT (1800) 

C0MM0H/BLK3/  MNO (84 1 ),NOD(84 1,4) 

COMMON /BLK5 /  SL ( 1 800 ) , SLP ( 1 800 ) , SLPP ( 1 800 ) 

C0MM0N/BLK8/  C(4,4 ),SO(4),XPV(5),YPV(5),XJCOB(5),C1 (8,8),ZY(8), 

•  FV(5,4),GV(5,4),XNV(5,4),SICT(4),DSIG(4),EPT(4),DEP(4), 

•  ST0R(6),PWPT,DPWP,CSC(4,4) 

DIMENSION  U(2),UX(2),UY(2) 

C 

C  FORMAT  STATEMENTS 
C 

920  F0RMAT(I8,3X,1P9E12.3) 

922  P0RMAT(//,5X,7HELEMENT,45X, 

1  28HELEMENT  STRAINS  AND  STRESSES, 6X, 3HN0. ,5Z| 

2  9HEPSIL0N-X, 3X,9HEPSIL0N-Y, 3X , 9HEPSIL0N-Z 

3  , 3X, 8HGAMI'A-xy,4X, 7HSIGMA-X, 5X, 7HSIGMA-Y, 5X, 7HSZCMA>Z , 

4  5X,6HTA0-XY,10X,’U’ ) 

923  F0RMAT(//,5X,’ELEMENT',45X, 

1  ’ ELEMENT  STRAINS  AND  STRESSES ' , / , 6X , ' NO . * , 5X , 

2  • EPSILON-R’ , 3X, ’ EPSILON -Z ' , 3X, ’EPS-THETA ' 

3  ,3X, 'GA»ltA-RZ',4X,’SlGMA-R',5X,»SICMA-Z’,4X,’SIG-THETA', 

4  4X, 'TAD-RZ',10X,'U') 

924  FORMAT (1H0,4X,4HN0DE,  7X,  13HDISPLACDJENTS,  /  6X,  3HW0.  ,  8X, 

1  1HU,  10X,  1HV) 

C 

C  FOR  EACH  ELEttENT  FIND  STRAINS  AND  STRESSES 
C 

IFdTSTOP  .EQ.  0)  GO  TO  110 
IF(MTYPE  .LT.  2)  HRITE(6,922) 

IF(KTYPE  .BQ.  2)  HRITE(6,923) 

110  DO  760  1X8 1, NELEM 
MN8MN0(IX) 

N0D(IX,1)=IABS(N0D(IX,1)) 

C 

C  RECALL  ELEMENT  INFORMATION  FROM  DISK 
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onno  o  o  ooo 


READ(2=IX)((S0(J),(C(I,J),CSC(I,J),I=1,H),Js1,4),(XPV(K),ypV(K), 

•  XJC0B(K),(FV(K,L),GV(K,L),XNV(K,L),Ls1,‘I),K=1,5),(SICT(M), 

•  DSIG(M),EPT(K),DEP(M),M=1,4),PMPT,DPWP,(ST0R(N),N=1,6)) 

CALCULATE  THE  STRESS  AKD  STRAIN  AT  THE  ELEMENT  CENTER 

DO  157  J=1,2 

DUUO.O 

DU2=0.0 

DU3=0.0 

DO  155  1=1, U 

NN=NOD(IX,I) 

II=NQ(NH)fJ-1 

••••••••••NOTE»»»»»«»NOTE«»»»«»NOTE 

UN=SL(II) 

DU1=DU1  +  XNV(5,I)»UN 
DU2=DU2  +  FV(5,I)«UN 
155  DU3=DU3  GV(5,I)«UN 
U(J)=DU1 
UX{J)=DU2 
157  Uy(J)=DU3 
R0=0.0 

IF(MTyPE  .EQ.  2)  RO=1.0/XPV(5) 

DEP(1)=UX(1) 

DEP(2)=Uy(2) 

DEP{>»)zux(2)  +  uy(i) 

DEP ( 3 ) =R0  *0 ( 1 ) -XKP« (CSC ( 1 , 3 ) •DEP ( 1 )+CSC (2,3) •DEP ( 2 ) 

1  +  CSC(3,*‘)«DEP(A)+  SO(3))/CSC(3,3) 

DO  ‘•00  1  =  1, A 
DU=SO(I) 

IF(NONLIN  .EQ.  2)DU=0.0 
DO  380  J=1,4 

380  DU=DU  +  CSC(I,J)»DEP(J) 

A 00  DSIG(I)=DU 

DSIG(3)=DSIG(3)«(1.0  -  XKP) 

DPVPrO.O 

IFdTSTOP  .EQ.  0)  GO  TO  750 

IF  COIIVERCENCE  HAS  OCCURRED,  SUM  ELEMENT  STRESSES  AND  STRAINS 
AND  PRINT  RESULTS 

DO  7A0  J=1,A 
EPT(J)=EPT(J)  +  DEP(J) 

7A0  SIGT(J)=SICT(J)  ♦  DSIG{J) 

PWPT=PWPT  ♦  DPWP 

WRITE(6,920)  IX, (EPT(J),J=1,A), (SICT(J),J=1,A),PWPT 
750  CONTINUE 

WRITE(2=IX)((SO(J),(C(I,J),CSC(I,J),I»1,A),J=1,A),(XPV(K),ypV(K), 

•  XJC0B(K),(FV(K,L),CV(K,L),XKV(K,L),L=1,A),Ks1,5),(SI0T(K), 

•  DSIG(K),EPT(M),DEP(M),M=1,4),PWPT,DPWP,(ST0R(N),N=1,6)) 
760  CONTINUE 

IFdTSTOP  .EQ.  0)  RETURN 


ooo  ooo  ooooo  oon 


DZSPLACEHENTS  SOtMED  AND  PRINTED 

VRITE(6,924} 

DO  785  Isl.NPT 
JJsNQ(I) 

IP(JJ  .EQ.  NQ(  1+1) )G0  to  765 
DISPLT(JJ)OISPLT(JJ)  +  SL(JJ) 
DISPLT(JJ+1)=DISPLT(JJ+1)  +  SL(JJ+1) 
WRITE (6 , 920 )  I ,DISPLT( J J ) .DISPLT ( JJ+1 ) 
785  CONTINUE 
RETURN 
END 


SUBROUTINE  STIFTJS (IX,LSTND,LT,NCOL,S ) 

THIS  SUBROUTINE  FORKS  THE  ELEMENT  MATRIX  FOR  A  QUADRILATERAL 
ELEMENT 

COMl ;ON /BLKO/  NPT . NELEM, NBPTC , MTYPE , IHISBF. ITKAX , RELAX , ERMAX , 

•  VITFAC .MODNEV. ITFAC , NONLIN , NVA Y. IREPET, ITNO, NSTP 

•  ,ITIM,TIMF,TIMB,XKP,XXX 
COMMON/BLK1/  PROP(3,21),FXA(3),FYA(3),ITYPA{3) 

COrmJ /BLK2/  X (900  ) ,  Y (900) ,  NQ  (90 1 ) ,  DISPLT (1 800  ) 

COMMON /BLK3/  MN0(8‘lt ),HOD(8»«1,A) 

COM»;OH  /BLKlI  /  NODB  (1 20  ) ,  BIV  ( 1 20 , 3 ) 

COMMON /BLK5 /  SL ( 1 8  00 } , SLP ( 1 6  00 ) , SLPP (1800) 

COHMON/BLK8/  C(4,i| ),SO(ll ),XPV(5), YPV(5),XJCOB(5),C1  (8,8),ZY(8), 

•  FV(5,«),GV(5,‘I),XN?(5,'»),SICT(»»),DSIC(A),EPT(1»),DEP(A), 

•  ST0R(6),PWPT,DPWP,CSC(H,A) 

DIMENSION  IIFLG(2).BIVD(2),S(LT,NC0L)»C2(8,6) 

RECALL  ELEMENT  INFORMATION  FROM  DISK 

READ(2=IX)((S0(J),(C(I,J),CSC(I,J),I»1,4),J=1,‘I),(XPV(K),YPV(K), 

•  XJC0E(K),(FV(K,L),CV(K,L),XNV(K,L),L=1,A),K=1,5),(SICT(M), 

•  DSIG(M),EPT(M),DEP(M),M=1,A),PWPT,DPWP,(ST0R(H),N=1,6)) 

CALCULATE  INCREMENTAL  PROPERTIES 

MNrMfJOdX) 

IN=2 

CALL  PROPTY(IN,MN,NZ,NZ) 

C11sC(1,1) 

C12=C(1,2) 

C13=C(1,3) 

CH»=C(1,4) 

C22=C{2,2) 

C23sC(2,3) 

C2iJ=C(2,»l) 

C33=C(3,3) 

C3'I=C(3,1») 

ClHJ»C(A,4) 

ClIXsCSCd,  1) 
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ooo  non  ooo  ooo  ooo 


C12X=CSC(1,2) 

C13X=CSC(1,3) 

CUX=CSC(1,4) 

C22XsCSC(2,2) 

C23X=CSC(2,3) 

C24X=CSC(2,4) 

C33XsCSC(3,3) 

C34XsCSC(3,4) 

C44X=CSC(4,4) 

SUSO(I) 

S2sS0(2) 

S3=SO(3) 

S4=S0(4) 


INITIALIZE  ELEMENT  MATRICES 


00  60  J=1,8 
DO  50  K:1,8 
C2(J,K)=0.0 
50  C1(J,K)=0.0 
60  ZY(J)=0.0 

SET  PARAMETERS  FOR  TYPE  OF  2-D  ANALYSIS  TO  BE  PERFORMED 

Rsl.O 

ROsO.O 

PLANE  STRESS  TERliS 

D11=XKP»C13»C13/C33 
D12=XKP»C13»C23/C33 
D14=XKP»C13*C34/C33 
D22=XKP»C23«C23/C33 
D24=XKP»C23«C34/C33 
D44=XKP»C34»C34/C33 
CALL  I NTP ( IH ISBF , TIMB , TIKF , ITIM , DF ) 

MNsKNOdX) 
n=FXAa3I)*DF 
FY=FYA (MN)»DF 

NUMERICAL  INTEGRATION  LOOP 
DO  240  N=1,4 

IF(KTYPE  .LT.  2)  GO  TO  152 
R=XPV(K) 

R0=1.0/R 


CALCULATE  ELE!IENT  MATRICES 

152  D0=XJC0E(N)»R 
DO  220  1=1,4 
I1=2»I-1 
I2slU1 
FI=FV(N,1) 
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GIsGV(N,I) 

XNIsXNV(N,Z) 

ELEMENT  LOAD  MATRIX 

ZT(I1)=ZY(I1)  ♦  DU»{XNl«(FX-RO»S3)  -  FI«(S1-IXP«C13/C33»S3) 

1  -  GI»(Sll  -  XKP»C3R/C33»S3)) 

ZY(I2)=ZY(I2)  +  DU«(XNI»FY  -  FI»(SA-IKP»C3‘»/C33*S3) 

1  -  GI«(S2  -  IKP»C23/C33»S3)) 

IPCMODNEW  .EQ.  1)  GO  TO  220 

DO  210 

J1=2»J-1 

J2sJU1 

FJ=FV(N,J) 

GJ=GV(N.J) 

XNJ=XNV(N.J} 

ELEMENT  STIFFNESS  MATRIX 

C2(I1,J1)=C2(I1,J1)  ♦  DO*{(C11X)»FI«FJ  ♦  C33X»RO»RO»XNI»XNJ 

1  ♦  (CAAX)*GI*GJ  ♦  C13X»RO*(FI«XNJ  ♦  FJ*XNI) 

2  ♦  (C1AX)«(GI«FJ  ♦  GJ«FI)  ♦  C3‘»X»RO«(XNI«GJ+XNJ«GI)) 
C2(I1,J2)=C2(I1,J2)  ♦  D0»((CAAX)»0I»FJ  ♦  (C12X)»FI»GJ 

1  ♦  (C1AX)«FI«FJ  ♦  C23X»RO»XNI«GJ  ♦  (C2AX)«OI«CJ 

2  ♦  C3<X«RO»XNI»FJ) 

C2(I2,J1)«C2(I2,J1)  ♦  DO»((CAAX)«FI»GJ  ♦  (C12X)»FJ«»GI 
^  *  (C1AX)*FI*FJ  ♦  C23X*R0»0I»XHJ  ♦  (C2AX)»0I»GJ 

2  ♦  C3‘»X»RO«XNJ»FI) 

C2(I2,J2)aC2(I2,J2)  +  D0»((C22X)»0I»CJ  ♦  (CA«IX)»n«FJ 
1  ♦  (CZAD^CFItKJJ  ♦  FJ«GI)) 

C1(I1,J1)=C1(I1,J1)  ♦  DU»((Cn-Dl1)»FI«FJ  ♦  C33»RO«RO»XNI«XNJ 

1  ♦  (C4A-D«I4)»GI«GJ  +  C13»RO»(FI»XNJ  ♦  FJ»XNI) 

2  +  (C14-D14)»(0I»FJ  ♦  GJ»FI)  ♦  C34»R0»(XNI»0J+XNJ«GI)) 
C1(I1,J2)=C1(I1,J2)  ♦  D0»((C1|I»-D44)»GI»FJ  ♦  (C12-D12)»FI»0J 

1  ♦  (C14-D14)»FI»FJ  ♦  C23«R0»XNI«»CJ  ♦  (C24-D24)«CI»GJ 

2  ♦  C34«R0»XNI»FJ) 

C1(I2,J1)=C1(I2,J1)  ♦  D0»((C44-D44)»FI»0J  ♦  (C12-D12)»FJ»GI 

1  ♦  (C14-D14)»FI»FJ  ♦  C23«RO»OI»XNJ  ♦  (C24-D24)«OI»GJ 

2  +  C34«R0»XNJ*FI) 

210  C1{I2,J2)=C1(I2,J2)  ♦  D0»((C22-D22)»GI»GJ  +  (C44-D44 )«FI»FJ 
1  +  (C24-D24)»(FI»GJ  ♦  FJ»GI)) 

220  CONTINUE 
240  CONTINUE 

DO  245  J=1,8 
DO  245  1=1,8 
C2(I,J)sC2(J,I) 

245  C1{I,J)sC1(J,I) 

IFCNONLIN  .EQ.  1 )GO  TO  270 

DO  260  1=1,8 

DUbO.O 

JaO 

DO  250  JJ=1,4 
DO  250  N=1,2 
L=NOD(IX,JJ) 


o  o  o  o  o  r> 


L=NQ(L)+N-1 

J=J+1 

250  DU=DU+C2(I,J)»SLP(L) 

260  ZY(I)=-DU 
270  CONTINUE 
C 

C  THE  NODE  POINT  SPECIFICATIONS  ARE  CONSIDERED 

C 

NU=2 

DO  33H  J=1,4 
NRQ=NOD(IX,J) 

KLsNBPTC  +  1 
DO  331  K=1,KL 
NR=NU*(J-1) 

IF(K  .LT.  KL)  GO  TO  300 

APPLY  BOUNDARY  CONDITION  WHEN  R=0 

IF(KTYPE  .LT.  2  .OR.  X(NRQ)  .NE.  0.0)  GO  TO  331 
IIFLG(1)=1 
IIFLG(2)=0 
BIVD(1)=0 
BIVD(2)=0 
GO  TO  320 

300  KK=IABS(HODB(K)) 

KIsKK/ 1000000 
IF  (NRQ  .NE.  K1)  GO  TO  331 
XX=1.0 

IF  (NODB(K)  ,LT.  0)  XXsO.O 
NODB(K)=-KK 

IIFLG( 1 )=MOD(KK, 100)/10 
IIFLG(2)=M0D(KK, 10) 

IH1=M0D(KK, 1000000)/10000 
CALL  INTP(IH1,TIMB,TIMF,ITIM,DF) 

BIVD(1)=BIV(K,1)»DF 
IH2=H0D{KK, 10000)7100 
CALL  INTP(IH2,TIMB,TIMF,IT1M,DF) 

BIVD(2)=BIV(K,2)-*DF 
AKC=BIV(K,3) 

IF(ANG  .EQ.  0.0)  GO  TO  320 

TRANSFORMATION  TO  LOCAL  COORDINATE  AXES 

CA=COS(ANG) 

SA=SIN(ANG) 

D1=ZY(NR.^1) 

D2=ZY(NR.»2) 

ZY(NR-f1)=D1»CA+D2«SA 
ZY  (NR.^2 ) = -D 1  •SA-^D2«CA 
DO  315  JJ=1,8 
D1=C1(MR*1,JJ) 

D2sC1(NR>2,JJ) 

C  1(IIR+ 1 ,  J  J )  r  D 1  •CA^D2«SA 
315  C1(NRf2,JJ)=-D1«SAfD2»CA 


D1=C1(1I,»R+1) 

D2sC1(II,NR4.2) 

C1(II,NR4-1  }8D1*CA4D2*SA 
318  C1(II,NR+2)=-D1»SA+D2*CA 
320  DO  330  Ns1.2 
DU1=BZVD(N} 

NRsKRt^l 

C 

C  LOADS  ADDED  IN 

C 

ZY(NR)=ZY(NR)>XX»Dai 
IF(IIFLG(N)  .BQ.  0)  GO  TO  330 
C 

C  DISPLACEMENTS  SPECIFIED 

C 

C«*«*««««cc«»*ff«»*||OTE***«*««**«*««*«*NOTE 

IFCNONLIN  .EQ.  1)G0  TO  325 

NRM=NQ(NRQ)-1i-N 

DOUDOI-SLP(NRM) 

325  ZY(HR)=DU1  •  XX  •  XXX 
C1(NR,NR)=XXX  •  XX 

330  CONTINUE 

331  CONTINUE 
33'»  CONTINUE 

C 

C  THE  ELEMENT  MATRIX  IS  NOW  ADDED  INTO  THE  SYSTEM  MATRIX 

C 

NRCCsO 
DO  355  K=1,l| 

KK=HOD(IX,K) 

NR=NQ(RK}-1 
DO  350  Mel, NO 
NRCCeNRCC-t-l 
NReNR  -1^  1 

IFCMODNEW  .EQ.  1)  GO  TO  350 
NRMeNR  -  LSTND 

NCCC=0 
DO  3'«5  L=1,A 
JJsNODdX.L) 

NCNeNQ(JJ)  -  NR 
DO  3*15  N=1,NU 
NCCC=NCCC+1 
NCN=NCN-*-1 

IF  (NCN.LT.  1)  GO  TO 
S(NRM,NCN)eS(NRM,NCN)  ♦  Cl (NRCC.NCCC) 

3A4  CONTINUE 
3A5  CONTINUE 

350  SL(NR)eSL(NR)  ZY(NRCC} 

355  CONTINUE 

HRITE(2=IX)((S0(J),(C(I,J),CSC(I,J),I=1,A),J=1,A),(XPV(K),YPV(K), 

•  XJCOB(K),(FV(K,L),OV(K,L),XMV(r,L),L=1,4),Ke1,5),(SIGT(M), 

•  DSIG(M},EPT(K),DEP(M),Me1,A),PVPT,DPHP,(ST0R(N),Ne1,6)) 


ooon  ooo  noooooo 


RETURN 

END 

SUBROUTINE  SOLVE  (L1,L2,LT,NR0W,NC0L,ID1SK,Q7,S) 

THIS  SUBROUTINE  FORMS  AND  SOLVES  THE  SIMULTANEOUS  EQUATIONS. 
VARIABLE  DIMENSIONING  IS  USED  TO  MAXIMIZE  THE  LENGTH  OF  THE  MAIN 
BLOCK.  DISK  STORAGE  IS  USED  WHEN  MORE  THAN  ONE  BLOCK  IS  REQUIRED. 
WHEN  M0DNEW=1,  ONLY  REDUCTION  OF  THE  RIGHTHAKD  SIDE  AND  BACK 
SUBSTITUTION  TAKES  PLACE. 

C0»1I  ION  /BLKO/  NPT ,  NELEM ,  NBPTC ,  MTYPE ,  IHISBF ,  ITMAX ,  RELAX ,  ERMAX , 

•  VITFAC , MODNEW, ITFAC , NONLIN , NWAY, IREPET, ITNO , NSTP 

•  ,ITIK,TIMF,TIMB,XKP,XXX 
COMMON /BLK2/  X (900 ) , Y(900 ) , NQ (90 1 ) , DISPLT ( 1600 ) 

COMMOH/BLK3/  MNO(84l ),MOO(81i  1  ) 

COMMON /BLR5 /  SL ( 1 8  00 ) , SLP (1 800 ) , SLPP ( 1 800 ) 

DIMENSION  S(LT,NCOL),Q7(1) 

THE  STIFFNESS  MATRIX  IS  GENERATED  IN  BLOCKS  AND  STORED  ON  DISK 
ID=0 

288  IDsID  *  1 

LSTND=(ID-1)»L1 
!IWND=LSTND  +  LI 

EACH  ELEMENT  IS  EXAMINED  TO  DETERMINE  IF  IT  CONTRIBUTES 
TO  THE  BLOCK 

DO  355  I=t, NELEM 
IX=I 

IF(NOD(I,1)  .LT.  0)  GO  TO  355 
DO  300  L=1,A 
KK=NOD(I,L) 

K1=NQ(KK) 

IF  (  K1  .LE.  NWND)  GO  TO  305 
300  CONTINUE 
GO  TO  355 

:  CALCULATE  THAT  PORTION  OF  THE  STIFFNESS  MATRIX  GIVEN  BY  A 

:  CONSIDERATION  OF  ELEMOJT  I 

305  CALL  STIFNS  (IX,LSTND,LT,MCOL,Q7) 

HOD(I,1)=-NOD(I,1) 

355  CONTINUE 

:  THE  BLOCK  OF  EQUATIONS  IS  REDUCED  AND  PUT  ON  DISK  IF  REQUIRED 

IF(MODNEW  .EQ.  0  )  GO  TO  400 

IFdDISK  .NE.  0)  READdsID)  ( (S(N,M),M»1,NCOL),Hs1,Ll ) 
or  TO  505 


*>*•'  JCE  .HE  LEFTHAHD  SIDE 


400  DO  500  Nsl.LI 
DIAG=S(N,1) 

IF(DIAG  .EQ.  0.0)  GO  TO  500 
IsH 

DO  475  Ls2,NCOL 

CsS(N,L)/DIAG 

IsI+1 

IF(C  .EQ.  0.0)  GO  TO  475 
JsO 

DO  450  K=L,NCOL 
J=J+1 

450  S(I,J)sS(I,J)  -  C»S(N,K) 

475  CONTINUE 
500  COTITINUE 


C*««»«»pOR  LINEAR  SYSTEMS  INCLUDE  FOLLOWING  STATEMENT 
C  IF (ID  .GE.  IDISK)GO  TO  505 

IFdDISK  .GT.  0)  VRITE(UID)  ((S(N,M),Hs1,NCOL),Ns1,L1) 
C 

C  REDUCE  THE  RIGHTHAND  SIDE 

C 


505  DO  520  N=1,L1 
DIAG=S(N,1) 

IFCDIAG  .EQ.  0.0)  GO  TO  520 
NRsN  *  LSTND 
D=SL(NR)/DIAG 
IsNR 

DO  510  L=2,NCOL 
IsI+1 

CfctccvifiF  CANNOT  OVER  FLOW  SL  IN  COM1!OH,  INCLUDE  NEXT  STATaOT 
IF (I  .GT.  NR0W)G0  TO  520 
510  SL(I)=SL(I)  -  S(N,L)»D 
520  CONTINUE 

IF(ID  .GE.  IDISK  .OR.  HODNEW  .EQ.  1)  GO  TO  600 
C 

C  THE  NUMBER  TWO  BLOCK  OF  EQUATIONS  IS  SHIFTED  INTO  THE 

C  NUI^iBER  ONE  POSITION 

C 

NsLI 

DO  530  Is1,L2 
N=N+1 

DO  530  J=1,HCOL 
S(I,J)=S(N,J) 

530  S(N,J)=0.0 

IF(L2  .GE.  LI)  GO  TO  600 
LL=L2  +  1 
DO  540  IsLL.LI 
DO  540  J=1,NCOL 
540  S(I,J)eO.O 

600  IF(ID  .LT.  IDISK)  GO  TO  286 
C 

C  BACK  SUBSTITUTION 

C 

IDsIDISK  *  1 
NRsIDISK«LU1 


o  o  o 


IPCIDISK  .EQ.  0)NR:NR0W-i-1 
610  1D=ID  -  1 

IF(ID  .LT.  IDISK)  READ(1=ID)  {(S(N,M),M=1,NC0L>,Ns1,L1) 

N=L1  1 

DO  750  M=1,L1 

HRzNR  -  1 

N=K  -1 

IP (NR  .GT.  NROW)GO  TO  750 
DUIzSL(NR) 

IP(S(N,1)  .EQ.  0.0)  GO  TO  750 
L=NR 

DO  725  K=2,NC0L 
L=L  1 

CANNOT  OVER  FLOW  SL  IN  COMMON,  INCLUDE  NEXT  STATEMENT 
IF(L  .GT.  NROW)GO  TO  727 
725  DU1=DU1-S(N,K)»SL(L) 

727  SL(NR)=DU1/S(N,1) 

750  CONTINUE 

IF(ID  .GT.  1)  GO  TO  610 

RETURN 

END 


SUBROUTINE  INTP  (I,T2,T1,ITIM,DF) 

SUBROUTINE  TO  INTERPOLATE  HISTORY  FUNCTION 

C0KM0H/BLK7/  FUN ( 10, 3) t FUNT ( 10, 3) ,NPTS (3,3) 

DIMENSION  TT(2),F(2) 

DFsI.O 

IF(I  .LT,  0)  RETURN 
IF(I  .GT.  0)  GO  TO  20 
IFdTIH  .GT.  1)  DFzO.O 
RETURN 

20  DF=FUHT(1,I) 

IPdTIM  .EQ.  NPTSd.D)  RETURN 

NPTS(I,1)=ITIM 

TT(2)=T2 

TT(1)=T1 

NP=NPTS(I,2) 

N=NPTS(I,3) 

FUHT(1,I)=0.0 

IF(TT(2)  .LT.  FUNT(NP,I))  NP=1 

DO  300  LL=1,2 

L=3-LL 

T=TT(L) 

DO  100  J=NP,N 
K=J 

IF(T  .LE.  FUNT (K+ 1,1))  GO  TO  200 
100  CONTINUE 

200  F(L)=FUN(K,I)  +  (FUH(K+1,I)-FUN(K,I))«(T-FUNT(K,I))/ 
•  (FUNT(K+1,I)  -  F0NT(K,I)) 

300  CONTINUE 
NPTS(I,2)=K 
DFrF(l)  -  F(2) 


FUNT(1,I)=DF 

RETURN 

END 

SUBROUTINE  GEOM ( IX , MN , ISIGN, SCOEF, SIGH , SIGV , U , XC . TC ) 

c 

C  THIS  SUBROUTINE  INITIALIZES  THE  STRESS  AND  PROPERTY  ARRAYS  AND 
C  CALCULATES  THE  ISOPARAMETRIC  TRANSFORMATION  AND  STORES  ON  DISK 
C 

C0M>!0H/BLK2/  X (900 ) , Y(900 ) ,NQ (901 ) , DISPLT(  1800 ) 

COMMON/BLK3/  MN0(641  ),NOD(8il1,4) 

C(MI0N/BLK6/  R0A(<I).SCA(4),ETA(A) 

COMMOH/BLE8/  C(A,A ),SO(A),XPV(5),ypV(5),XJCOB(5),C1 (8,8),ZY(8), 

•  FV(5,‘»),GV(5,A),XNV(5,9),SIGT(ll),DSIC(ll>,EPT(«l),DEP(l»), 

•  STOR(6),PWPT,DPWPtCSC(4,ll> 

DIMENSION  SC0EF(8,6},RQ(4),ZQ(4) 

C 

C  THE  AND  NODE  POINT  COORDINATES  ARE  FOUND 

C 

DO  100  J=1,l» 

K=NOD(IX,J) 

RQ(J)=X(K) 

100  ZQ(J)=Y(K) 

C 

C  ISOPARAMETRIC  TRANSFORMATION  FACTORS 

C 

A1=ZQ(1)+ZQ(2)-ZQ(3)-ZQ(4) 

A2sZQ(1)-ZQ(2)-ZQ(3)-^ZQ(‘») 

B=  ZQ(1)-ZQ(2)+ZQ(3)-ZQ(*») 

C1P=RQ(1)-RQ(2)-RQ(3)+RQ(M) 

C2=RQ(1)*RQ(2)-RQ(3)-RQ(*) 

D=  RQ(1)-RQ(2)^RQ(3)-RQ(‘») 

C 

C  QUADRATURE  POINT  LOOP 

C 

DO  155  N=1,5 

IF(N  .LT.  5)  GO  TO  111 

SC=0.0 

BTrO.O 

GO  TO  IIH 

111  D0s1.0/SQRT(3.0) 

SC=DU»SCA(N) 

ET=DU*ETA(N) 

1U»  D1=C1P+D»ET 
D2rAUB»SC 
D3=C2+D*SC 
DA=A2+B»ET 
DU=1.0/(D1»D2-D3»DA) 

C 

C  CALCULATION  OF  SHAPE  FUNCTION  DERIVATIVES 

C 

XCrO.O 

YC=0.0 

DO  150  Is1,4 
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DS^ROAd) 

D6=SCA(I) 

D7=ETA(I) 

FV(N,I)*DU«(D2»(D6  ♦  D5»ET)  -  DM«(D7  ♦  D5«SC>) 

GV(N,I)=DU«(D1*(D7  ♦  D5«SC)  -  D3»(D6  D5»ET)) 

D8=0.25»(1.0  +  D6«SC)«(1.0  ♦  D7»ET) 

XNV(N,I)=D8 
XCrXC  +  HQ(I)»D8 
150  YC=YC  +  ZQ(I)*D8 
XPV(N)=XC 
YPV(N)=YC 
XJC0B(N)=0. 0625/DC 
155  CONTINUE 

INITIALIZE  STRESSES  AND  STRAINS 

0=0.0 
SIGHsO.O 
SIGV=0.0 
DPWP=0.0 
DO  160  1=1, A 
DSIG(I)=0.0 
DEP(I)=0,0 
160  EPT(I)=0.0 

IFdSIGN  .EQ.  0)  GO  TO  162 
SIGVsSCOEFdSIGN.I)  ♦  SC0EFdSIGH,2)»YC 
SIGH=SC0EFdSIGH,3>  +  SC0EFdSI0H,A)*YC 
UsSC0EFdSIGN,5)  ♦  SC0BFdSIGR,6)«YC 
162  SIGTCDsSIGH 
SIGT(2)=SI0V 
SIGT(3)=SIGH 
PWPT=U 
SIGT(A)=0.0 

STORE  INFORMATION  OH  DISK 

WRITE(2=IX)((S0(J),(Cd,J),CSCd,J>,I=1,A),J=1,A),(XPV(K),YPV(K), 

•  XJCOB(K),  (FV(K,L),GV(K,L),XNV(K,L),L=1,A),Ks1,5),  (SICT(M), 

•  DSIG{K),EPT{M),DEP(M),Ms1,‘J),PWPT,DPWP,(ST0R(H),Ns1,6)) 
RETURN 

END 


SUBROUTINE  RPROP(PROP) 

THIS  SUBROUTINE  READS  IN  AND  SCALES  THE  PROPERTIES  REQUIRED 
BY  THE  BOUNDING  SURFACE  PLASTICITY  MODEL  FOR  COHESIVE  SOILS. 

DII5ENSIOH  PROP (21) 

READ(5,801)  (PROPd),I=1,3),(PROPd),Is9,11),PROP(7), 

1  PR0P(21),PR0P(17),PR0P(16),PR0P(20),PR0P(5), 

2  PR0P(6),PR0P(8),PR0P(4),PR0P{18),PR0P(15), 

3  PR0P(12),PR0P(13),PR0P(19),PR0P(1A) 

801  F0RMAT(8E10.3) 

HRITE(6,901)  (PR0Pd),I»1,3).PR0P{A) 
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901  FORMAT (5Z, 'CLAY  PROPERTIES' //15Z, 'LAMBDA  r ' .E10.3/15X, 

1  'EAPPA  s',E10.3/15X,'MC  *' ,E10.3/15X, 

2  'ME/MC  s',E10.3) 

WRITE(6,902)  PROP{9),PROP(12),PROP(10),PROP(13), 

1  PROP(11),PROP(7),PROP{21),PROP(17), 

1  PR0P(19),PR0P(16),PR0P(18), 

1  PROP<HI),PROP(20) 

902  FORMAT (15X,'RC  =* ,E10. 3, 15X, 'RB/RC  «*,E10.3/ 

1  15X,'AC  =',E10.3,15X,'AE/AC  »',E10.3/ 

2  15X,'T  =',B10.3,15X,'PL  ■•,E10.3/ 

3  15X,'P0  =',E10.3/15X,'BARDENI)IG  PARAMETERS:'/ 

A  19X,'SMC  =',B10.3,10X, 'SME/SMC  s',E10.3/ 

5  19X,'HC  *',E10.3,10X,'HB/HC  s', £10.3/ 

6  ISX.'PROJECTIOM  POIMT  1/C  =',E10.3/ 

7  1 5X , ' INITIAL  VOID  RATIO  s  * , E 1 0. 3 ) 

IF(PR0P(5).LT.0.5)  WRITE(6,903)  PR0P(5) 

IF(PR0P(5).GE.0.5)  WRITB(6,90«)  PR0P(5) 

903  FORMAT(15X,'POISSOn«S  RATIO  »',E10.3) 

90A  FORMAT ( 1 5X , ' SHEAR  MODULUS  s ' , E 1 0. 3 ) 

WRITE (6, 905)  PROP (8) 

905  FORMAT (15X, 'ATMOSPHERIC  PRESSURE  s',E10.3) 

WRITE (6, 906)  PROP (15) 

906  FORMAT (15X, 'SIZE  OF  ELASTIC  ZONE  =',E10.3//) 
IF(PR0P(6).EQ.0.0)  WRITE(6,907) 

907  P0RMAT(5X, '•••••  DRAINED  CONDITIONS  •••••'//) 
IF(PR0P(6).NE.0.0)  WRITE(6,908)  PROP(6) 

908  FORMAT(5X, ••••*•  UNDRAINED  CONDITIONS  -  THE  COMBINED  ', 

1  'SKELETON  AND  HATER  BULK  MODULUS  s*,E10.3//) 

PROP(3)=PR0P(3)/(3.0«SQRT(3.0)) 

PROP(7)=PROP(7)«3.0 

PROP(21)=PROP(21)»3.0 

RETURN 

END 


SUBROUTINE  CLAY ( IDIM, INC, ITNO, PROP . STOR , SIGBH , EPM, 

1  DSIGM , DER] , C , CB , UB , DLTAU , GAM , KIND . LARGE ) 

C 

C  SUBROUTINE  TO  EVALUATE  YANNIS  DAFALIAS*  BOUNDING 
C  SURFACE  PLASTICITY  MODEL  FOR  CLAY  SOILS.  PREPARED  BY 
C  L.R.  HERRMANN  AT  THE  UNIVERSITY  OF  CALIFORNIA,  DAVIS  CAMPUS. 
C 

DIMENSION  PROP(21),STOR(6),SIGB(6),DSIG(6),DEP(6),C(6,6), 

1  SB(3,3),SF(3,3),II(6),DLTA(3,3),DEPM(6), 

2  SIGBK(6),DSIGM(6),DEPT(3,3)>EPH(6),EPB(6),CB(6,6) 
DATA  11/11,22,33,12,13.23/,  DLTA/1.0, 3*0.0, 1.0, 3*0.0, 1.0/ 
ALFUN (CV , RT , SINV ) =2. 0»RT*CV/ ( 1 . 0+RT- ( 1 . 0-RT ) •SINV ) 
SMALL=0.0001*PROP(8) 

DO  40  Is1,6 
SIGB(I)sSIGBM(I) 

DSIG(I)sDSIGM(I) 

EPBCI)=EP1!(I) 

40  DRP(I)xDBPM(I) 

IF(ITNO.OT.I)  GO  TO  100 
IFCINC  .0T.1)  GO  TO  50 


ooooooo  ooo  ooo  ooo 


INITIALIZE  HISTORY 

ST0R(1)spR0P(21) 

ST0R(2)=ST0R(1) 

STOR ( 3 ) = 0 . 5» (SIGB ( 1 )+SIGB  <  2 ) ) 

STOR('»)=0.01»PROP(8) 

STOR(5)=0.0 
GO  TO  100 

UPDATE  HISTORY 

50  ST0R(1)=ST0R(2) 

STOR ( 3 )sST0R ( 3 )+ST0R (« ) 

STOR (5 ) sSTOR (5 )^ST0R (6 ) 

CONVERT  FROM  PLANE  STRAIN  TO  3-D 

100  IF(IDIM.EQ.3)  GO  TO  200 
SIGB(4)=SIGB(3) 

SIGB(3)=STOR(3) 

DSIG(A)=DSIG(3) 

DSIG(3)=ST0R(1») 

DEP(A)=DEP(3) 

DEP(3)s0.0 

EPB(l|)rEPB{3) 

EPB(3)sO.O 
DO  110  Is5,6 
SIGB(I)=0.0 
DSIG(I)s0.0 
EPB(I)s0,0 
110  DEP(I)=0.0 

DETERMINE  3-D  INCREMENTAL  PROPERTIES 
200  GAM=PR0P(6) 

CALCULATE  EFFECTIVE  STRESS  INVARIANTS  AND  DISTORTIONAL  STRESS 
AND  CHANGE  MATRIX  COMPONENTS  TO  TENSOR  COMPONENTS. 

XIB=0.0 

XIF=0.0 

DDILsO.O 

DILB=0.0 

DO  205  1=1,3 

DDIL=DDIL+DEP(I) 

DILB=DILB-*-EPE(I) 

XIB=XIB+SIGB(I) 

205  XIF=XIF<^SIGB(I)-»DSIG(I) 

VOIDB=1.0<^PROP(20) 

VOIDFrVOIDB 

IF (LARGE. EQ.O)  GO  TO  210 
VOIDB=VOIDB«EXP ( -DILB } 

VOIDFs V0IDF»EXP ( -DILB-DDIL ) 
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210  DO  215  Nc1.6 
I«II(N)/10 
JrMOD(II(N},10) 

SB(I,J)sSIGB(N)>ZIB«OLTA(I,J)/3.0 

SB(J,I)sSB(I,J) 

DBI>T(I.J)sDEP(N)*(1.0«DLTA(l,J))«0.5 

DBPT(J,I)=DEPT(I,J) 

SF  ( I ,  J }  sSIGB  (N  )-^DSlG  (N  )-DLTA  (I ,  J  )*ZIP/3. 0 
215  SF(J,I)=SF(I,J) 

GA)IP=0.0 

IF(K1ND  .EQ.  0)G0  TO  217 

GAMP=GAM 

0B=ST0R(5) 

DLTAUsGAM«DDIL 
217  XIBsXIB-UB»3.0 

XIFrXIF- (UB+DLTAO )«3.0 

ST0B(6)=DLTAU 

SRTJBsO.O 

SRTJFrO.O 

DO  220  1=1,3 

DO  220  J=1,3 

SRTJB=SRTJBfSB(I.J)«SB(I,J) 

220  SRTJF=SRTJF+SF(I,J)«SF(I,J> 

SRTJBsSQRT (0. 5«SRTJB) 

IF(SRTJB»1000.  .LT.  XIB)SRTJB=0.0 
SRTJFsSQRT (0. 5»SRTJF) 

IF(SRTJF*1000.  .LT.  1IF)SRTJF=0.0 

SCOBsO.O 

SCUFsO.O 

DO  225  1=1,3 

DO  225  J=1,3 

DO  225  K=1,3 

SCUB=SCUB+SB ( I , J ) tSB ( J , K )»SB (K , I ) 

225  SCOF=SCUF+SF ( I , J )»SF( J , K )«SF (K, I ) 

SCUB=SCUB/3.0 

SCUF=SCUF/3.0 

SN3AB=0.0 

IF (SRTJB . GT. SMALL )  SN3ABc 1 . 5*SQRT ( 3. 0 }*SCDD/SRTJB«»3 
IF(SN3AE.GT.  1.0)  SN3AB=  1.0 
IF{SN3AB.LT.-1.0)  SN3AB=-1.0 
SN3AF=0.0 

IF (SRTJF. OT . SMALL )  SN3AF» 1 . 5»SQRT (3.0) •SC0F/SRTJF»«3 
IF(SH3AF.GT.  1.0)  SN3AFs  1.0 
IF(SN3AF.LT.-1.0)  SN3AF=-1.0 
CS3AB=SQRT( 1 . 0-SN3AB*«2) 

CS3AF=SQRT ( 1 . 0-SN3AF»*2 ) 

AVOID  ZERO  MEAN  PRESSURE 

IF(ABS(XIB).GT. SMALL)  GO  TO  227 

DOxXIB 

XIB=SMALL 

IP(DU.LT.O.O)  XIB=>SMALL 
227  IF(ABS(XIF).GT.SMALL)  GO  TO  230 


ooo  non  non 


DU=X1F 

XIFrSMALL 

IF(DU.LT.O.O)  XlFs-SMALL 
230  CONTINUE 

CALCULATE  ELASTIC  PROPERTIES 

D01=VOIDB/3.0/PROP(2) 

DU2=1.5»(1.0-2.0*PROP(5))/(1.0+PROP(5)) 

DU=XIB 

IF(DU.LT.PR0P(7))  DU=PR0P(7) 

BB=DU1»D0 

GB=DU2»BB 

IF(PR0P(5).GT.0.5)  GB=PROP{5) 

DU  1 = VOIDF/3 . 0/PRO*'  C  2 ) 

DU=XIF 

IF(DU,LT.PR0P(7))  DU=PROP(7) 

BF=DU1«DU 

GF=DU2»BF 

IF(PR0P(5).GT.0.5)  GF=GB 
DO  235  M=1,6 
I=II(H)/10 
J=MOD(II(M), 10) 

DO  235  N=H,6 
K=1I(H)/10 
L=M0D(II(N), 10) 

DU  1 =DLT A ( K , I ) •DLTA ( L , J ) ♦DLTA ( K , J ) ‘DLIA ( 1 , L ) 

C(M, N)=GF*DUU(BF+G  AMP-2. 0*GF/3.0)»DLTA(I,J)»DLTA(K,L) 

CB(M,N)=GB«DUU(BB+GAKP-2.0«CB/3.0>«DLTA(I,J)«DLTA(K,L) 

CE(N,M)=CB(M,N) 

235  C(N,M)=C(F.,N) 


CALCULATE  SIZE  OF  BOUNDING  SURFACE 

XIOErSTORf 1) 

XI0FiST0R(2) 

XIL=PROP(7) 

DU10=1.0/(PR0P(1)-PROP{2)) 
IF(XIOB.GE.XIL.AND.XIOF.GE.XIL)  CO  TO  2«I0 
XIOBS=XIOB 

IFCXIOB.LT.XIL)  XIOBS=XIL 
XIOFS=XIOF 

IF(XIOF.LT.XIL)  XIOFSsXIL 

XIOF= XIOB+DU 1 0*0 . 5 • { ( XlOFS*VOIDF+XIOBS»VOIDB ) •DDIL- 
1  ( XIOBS»VOIDB/DB+XIOFS»VOIDF/BF ) • (XIF-XIB ) /3 . 0 ) 

GO  TO  2M5 

2M0  XIOF=XIOD»£XP(DU10«0.5«((VOIDB+VOIDF)»DDIL- 
1  (VOIDB/BB+VOIDF/BF)»(XIF-XIB)/3.0)) 

2U5  ST0R(2)=XI0F 

IF(INC+ITN0.EQ.2)  GO  TO  UlO 

CALCULATE  BOUNDING  SURFACE  PROPERTIES 

CALL  BOUND ( PROP , SRTJB, SN3AB, XSB , XIOB, XIB , QAHB, DFIB, 
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1  DPJB.X1CSB,DFALB,DFJJB,BSB,?0IDB) 

CALL  B0UND(PIK>P.SRTJF,SN3AF,XSF,XI0F,ZIF,GAMF,DP1P, 

1  DFJF.XKSF.DFALF.DFJJF.BSFtVOIDF) 

DBrBSB-1.0 

IFCDB  .LT.  0.0)  DBsO.O 
DF=BSF-1.0 
IF(DF.LT.O.O)  DFsO.O 

CALCULATE  PLASTIC  HOOOLUS 

CHECK  FOR  ELASTIC  ZOHE  AND  UNLOADING 

XMS= ALFUN ( PROP (17), PROP (19 ) , SN3AB ) 

DU7=0.0001M(1.0/XMS) 

LBrO 

DDD=  1 .  Oi-DD*  ( 1 . 0-PR0P(  15  )  } 

IF(DDD.LE.O.O)  GO  TO  352 
LB:1 

HsALFUN ( PROP ( 1 6 ) , PROP ( 1 8 ) , SN3AB ) 

DOsABS(XSB) 

IF(DU.LT.0U7)  DDsOU7 
D08=9. 0»DFIB»*2+DFJB««2/3. 0 
DU9sXI0B 

IF(XIOB.LT.XIL)  D09*XIL 

XKB=XKSB-fH«DB/DDD*  ( 1.  (W1. 0/D0*«XMS  )*D08«DU9*DU  10«V0IDB 
DO Is 3. 0»BB»DFIB 
DU2sGB«DFJJB 
DU2PsSQRT( 3. 0)«GB»DFALB 

D03sXKB+9.0»BB»DFIB»»24GD»DFJB»«2>»CB«(DPALB»CS3AD)»*2 

SUMsO.O 

TIsO.O 

IF(SRTJB«»2  .EQ.0.0)  go  to  350 

DO  3'»0  Is1,3 

DO  3*0  Js1,3 

DUsO.O 

DO  330  K=1,3 

330  DUsDUt^SB(I,K)«SB(K,J) 

T 1  =T  1  ♦  ( DO- 1 . 5»SCDB«SB  ( I ,  J ) /SRTJB”2 )  tDEPT  ( I ,  J  ) /SRTJB»»2 
3*0  SOMsSOM+SB(I,J)»DEPT(I,J) 

T1sT1-2.0*DDIL/3.0 

350  DUs  (D0 1  •DDIL*D02»S0M+D02pm  )/DU3 
IF(DO.LT.O.O)  LBsO 
352  LF=0 

DDDs 1 . 0+DF» ( 1 . O-PROP ( 15 ) ) 

IF(DDD.LE.O.O)  GO  TO  358 
LFsl 

HsALFUN(PROP(16),PROP(18),SH3AP) 

DUsABS(XSF) 

IF(D0.LT.D07)  D0sD07 
D08s9.0«DFIF»«2+DFJF»»2/3.0 
XKSsALFUN (PROP( 17 ) , PROP( 19 ) ,SN3AF) 

D09sXI0P 

IF(XIOF.LT.XIL)  DU9»XZL 

XKFsXXSF4^H«DF/DDD« ( 1 . 0f1 . 0/DU««XKS )*DU8«DU9«D0 10«VOIDF 
DU*s3.0«BF«DFIF 


D05=GF«DFJJF 

DU6=XKF+9.0«BF*DFIF**2*GF»DFJF»*2fOF»(DFALF»CS3*F)*»2 
D05P=GF»DFALF»SQRT ( 3 • 0 ) 

SOHsO.O 

11=0.0 

IF(SRTJF*«2  .EQ.0.0)  GO  TO  357 
DO  356  1=1,3 
DO  356  J=1,3 
DD=0.0 

DO  355  K=1,3 

355  DU=DUfSF(I,K)»SF(K,J) 

T 1 =T 1 ♦ (DO- 1 . 5»SCUF«SF (I , J ) /SRTJF»»2 ) •DEPT ( I , J ) /SRTJF»*2 

356  SUM=SUM+SF(I,J)«DEPT(I,J) 

T1=T1-2.0»DDIL/3.0 

357  DU=  (DU»«*DDIL+DU5»SUM+D05P»T1  )/D06 
IF(DU.LT.O.O)  LF=0 

CALCULATE  PLASTIC  PORTION  OP  INCREMENTAL  PROPERTIES 

358  IF(LF+LB.EQ.O)  GO  TO  A 10 
DO  400  M=1,6 
I=1I(M)/10 
J=M0D(1I(M),10) 

DO  400  N=M,6 

K=II(N)/10 

L»M0D{II(N),10) 

D0=0.0 

IF(LB.EQ.O)  GO  TO  373 

T2=0.0 

TUO.O 

IF(SRTJB»»4  .EQ.0.0)  CO  TO  370 
DO  360  LL=1,3 
T2=T2+SB(K,LL)»SB(LL,L) 

360  T1=TUSB(1,LL)»SB(LL,J) 

T 1 =DU2P* (T 1 /SRTJB»«2- 1 . 5»SCUB»SB ( I , J ) /SRTJB»»4-2 . 0»DLTA ( I , J ) / 3 . 0 ) 
T2=DU2P« (T2/SRTJB«»2- 1 . 5«SCUB»SE (K , L ) /SRTJB»«4-2. 0»DLTA (K , L ) /3 . 0 ) 
370  DUBs-(DU1»DLTA(I,J)+D02*SB(I,J)+T1)«(DU1»DLTA(K,L)+ 

1  DU2»SB(K,L)fT2)/D03 

IF(LF.EQ.O)  CO  TO  396 
373  T2=0.0 
71=0.0 

IF(SRTJF«»4  .EQ.0.0)  00  TO  390 

DO  380  LL=1,3 

T2=T2^SF ( K , LL ) ‘SF (LL , L ) 

380  T1=T1+SF(I,LL)«SF(LL,J) 

T1=DU5P»(T1/SRTJF«»2-1.5»SCUF»SF(I,J)/SRTJP«»4-2.0«DLTA(I,J)/3.0) 
T2sDU5P* (T2/SRTJF»*2-1 . 5»SCUF»SF (K, L ) /SRTJF»»4-2. 0»DLTA (E, L )/3. 0 ) 
390  D0=-(DU4»DLTA(I,J)^D05»SF{I,J)+T1)»(DU4«DLTA(K,L)+ 

1  D05»SF{K,L)+T2)/DU6 

396  C(M,N)=DU-*^C(K,N) 

CB(M,N)=DUB-^CB(M,N) 

CB(N,M)=CB(K,N} 

400  C(N,H)=C(M,N} 

410  CONTINUE 


ooo  ooooo  ooo 


IP(IDIM.EQ.3)  RETORN 

CONVERT  3-D  PROPERTIES  TO  PLANE  STRAIN 

DOsO.O 
DO  A20  1=1, A 
DD=C(3,I)*DEP(I)-*^D0 
C(3,I)=C(A,I) 

A20  C(A,I)sO.O 
DO  A30  1=1,3 
C(I,3)=C(I,A) 

A30  C(I,A)sO.O 
STOR(A}sDO 
RETURN 
END 


SUBROUTINE  BOUND (PROP,SRTJ ,SN3A , X, XIO,XI , GAM, DPI ,DFJ , 
1  XKS,DPAL,DPJJ.BS,VOID} 

SUBROUTINE  TO  EVALUATE  RELATIONSHIP  OP  STRESS  STATE 
TO  BOUNDING  SURFACE 

DIKENSIOn  PROP(21),FSS(3) 

ALFON  ( CV ,  RT ,  SINV ) =2. 0«RT*CV/  ( 1 .  O^-RT-  ( 1 . 0-RT  )  •SINV  ) 
DFOH ( FUN , RT , FUNC ) =FUN»»2« ( 1 . O-RT ) / ( 2. 0«RT»F0WC ) 
XN=ALFUN(PROP(3),PROP(A),SN3A) 

DNALsDFUN (XN , PROP ( A ) , PROP( 3 ) ) 

R=ALFUN ( PROP (9 ) , PROP ( 12 ) ,SN3A ) 

DRALsDFUK (R , PROP (12), PROP{9 ) ) 

A= ALFUN ( PROP (1 0 ) , PROP ( 1 3 ) , SN3A ) 
DAAL=DFUN(A,PR0P(13)»PR0P(10)) 

YS=R»A/XN 

CCsPROPdA) 

SHIFT  PROJECTION  POINT 
D1*XI-XIO»CC 

IF(ABS(D1}.LT.0.001}  DIsO.OOl 

D2sCC>1.0/R 

D3=D1«D2 

D5=CC»(CC-2.0/R) 

Q  =SRTJ/D1 
QC=XN/(1.0-R«CC) 

Q0=1.0E<»20 

IF(CC.HE.O.O)  Q0=XN»(SQRT(1.0+yS«yS)-(1.0+YS))/R/CC 
IF(SRTJ.NE.O.O)  GO  TO  3 
IF(DI.OT.O.O)  GO  TO  10 
GO  TO  30 

3  IF(CC.LT.1.0/R)  GO  TO  5 
IF(Q  .GE.0.0)  GO  TO  10 
IF(Q  .LE.  QC)  GO  TO  10 
IF(Q  .GE.  00}  GO  TO  30 
GO  TO  20 

5  IP(Q  .GE.  QC)  GO  TO  20 
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IF(Q  .GE.0.0)  GO  TO  10 
IF(Q  .LE.  00)  GO  TO  20 
GO  TO  30 
C 

C  PROJECTION  ON  ELLIPSE  1 

C 

10  DH=D1»DU((R-1.0)»SRTJ/W)»»2 

BS=XI0«  ( -D3+SQRT  (D3«D3-M»  (D5+  ( 2 . 0-R  )  /B ) ) )  /M 
LOST=1 
GO  TO  100 
C 

C  PROJECTION  ON  HYPERBOLA 

C 

20  D6rSRTJ»(1.0/R+A/XN)/XN 
D7=D3+D6 

D8=D1«D1-(SRTJ/XN)”2 
BS=-0. 5»XI0» (D5-2. 0»A/R/XH )/D7 
IF(D8.EQ.C.O)  GO  TO  100 

BS=XI0«(-D7+S0RT(D7«D7-D8»(D5-2. 0«A/R/XN ) ) )/D8 
L0ST=2 
GO  TO  100 
C 

C  PROJECTION  ON  ELLIPSE  2 
C 

30  TspROPdl) 

POPsXN/SQRT ( 1 . 0+YS««2 ) 

XJO* A» ( 1 . 0+YS-SQRT ( 1 . 0*YS«»2) ) /YS 
BT»T« (XJO-T«FOP ) / (XJO-2. 0«T»FOP ) 
ROs(BT-T)/POP/XJO 
PSI=1.0/(R»(ET-T)) 

D9=T-BT+CC 

D10aD1»D9 

D1 1=D1«Dl4R0»SRTJ«SRTJ 

BS=XIO« ( -D 1 04SQRT (D 1 0*D 1 0-D 1 1 • (D9»D9-BT»BT ) ) ) /D 1 1 
LOST* 3 

100  XIBARsBS*(XI-XIO«CC}4XIO*CC 
IF(XIBAR.EQ.O.O)  XIBAR< 1 . OE-20 
TH=BS«SRTJ/XIBAR 
XsTH/XN 
DUsXIO 

IF(XI0.LT.PR0P(7))  DUsPR0P(7) 

DUS= 1 2 . 0»V0ID/ ( PROP ( 1 )-PBOP { 2 ) ) •XI0»»2»D0 
GO  TO  (110, 200, 300), LOST 
C 

C  NORMAL  CONSOLIDATION  ZONE 
C 

110  PSI=YS/(R-1.0)»«2 

D0sR«  ( 1 . 0+Xn+R*  (R-2. 0)»X*X ) 

OAMs(1.04(B-1.0)»SQRT(1.04R«(R-2.0)«X»X))/DO 

DFI»2.0«XIO«(GAM-1.0/R)«PSI 

DFJJ«2. 0«XIO»CAM*( (R-1 . 0}/XN )««2*PSI«BS/XIBAR 

DFJ«DPJJ«SRTJ 

XKSsDUS* (GAM- 1 . 0/R ) • (GAM^R-E. 0 )«PSI«PSI /R 
DFALsPSI •e . 0» (R- 1 . 0 )«TH»CAM»XIO« ( ( (R- 1 . 0 ) / (R»»2» 
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1  (2. 0/R-GAK- 1 . 0  )  )^1 . 0  )«DRAL-  (R- 1 . 0  )»DNAL/XN )  /» ••R 

RBTORN 

OVERCONSOLIDATED  ZONE 
200  DUs1.0-Z*(1.04^IS} 

0AK=-  (DU-^SQRT  ( ( X-YS- 1 . 0  )«»2+  (I»X- 1 . 0  )«YS«YS ) )  /  ( R«  (X»X- 1.0)) 
DPIs2. 0*XIO« (GAM-1 .0/R) 

DPJs2. 0*X10*  ( ( 1 .  Oi^YS  )/R-X«GAM  )/XN 
XXSsDOS*  (GAM- 1 . 0/R  )•  (D0*GAHi^2. 0«A/XN  )/R 
DPJJsDFJ/SRTJ 

DPALs6. 0»XI0« (DNAL» (TH»GAM/XN-1 . 0/R+A/(R«TH»GAM )-2. 0*A/XN ) / 

1  XN»«2+DRAL«( 1.0/TH-1.0/XN+A/(XN»TH»GAM))/R»»2+DAAL»( 

2  1.0/XN-1.0/(TH«GAM«R))/XN) 

RETURN 

TENSION  ZONE 

300  GA^:s  ( -T+BT-SQRT (BT«BT-R0»TH«TH»T«  (T-2. 0»BT ) ) )/ ( 1 . 0+RO»THnH ) 
DPI=2. 0«PSI«XIO»(GAM+T-BT) 

DFJJs2. 0*PSI*XI0*GAM«R0*BS/XIBAR 
DFJsDFJJ^SRTJ 

XKSsDUS*PSI«PSI«  (GA^:■^T-BT  )•  (GAM*  (BT-T  )+T»  (2. 0»BT-T  ) ) 
DYSAL=YS«  (DRAL/Ri-DAAL/A-DKAL/XN ) 

DPOPALsFOP* (DHAL/XN-YS»DYSAL/ ( 1 . 0+YS»YS ) ) 

DJOALxXJO*  (DAAL/A-DYSAL/YS  )*>^A«  ( 1 . 0/YS-FOP/XH  )«DYSAL 
DBTALs ( (T-BT )«DJ0AL-(T-2. 0«BT)»T»DFOPAL)/(XJ0-2. 0*T*P0P) 
DROALxDBTAL/POP/XJO-RO*  (DFOPAL/FOPi-DJOAL/XJO  ) 

DFALs3. 0»PSI»XIO»TH»GAM*(DB0AL+2. 0*R0»DBTAL/ (T+GAM-2. 0»BT) ) 

RETURN 

END 


SUBROUTINE  ACCEL(X2,X1,X,C,XL) 

THIS  SUBROUTINE  CALCULATES  THE  AITKENS  CONVERGENCE  FACTOR 
C=1.0 

DUs-X2+2.0»X1-X 

IF(DU  .EQ.  0.0)RETURN 

C=(X1-X2)/D0 

IF(C  .GT.  XDCxXL 

IF(C  .LT.  1. 0/XL )Cx 1.0/XL 

RETURN 

END 


/ 
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